Алгоритм поиска основной частоты, Метод максимального правдоподобия Опции

[QuickNick]

Здравствуйте, товарищи.

Мне нужно улучшить алгоритм поиска основной частоты по методу максимального правдоподобия. Был пример у меня, как делать, я его модифицировал в нескольких местах (добавил настраиваемость оконной функции и фильтра). Возможно, некоторые места вам покажутся совершенно абсурдными - но я пока новичок в DSP, некоторые места в старом алгоритме вызывают у меня сильное удивление, но пока не хватает опыта для однозначного решения: или я не понимаю момент, или предыдущий программист допустил косяк.

- Вход: массив байтов byteData длиной sampleLength, количество разбиений sampleRate.
- Выход: фундаментальная частота fundamentalFrequency и громкость fundamentalVolume.

Алгоритм:
1. Массив байтов преобразуется в вешественный массив амплитуд.
2. Проверяется максимальная амплитуда - если она ниже определённого значения, то это квалифицируется как почти полная тишина (выход из алгоритма)
3. Получение спектра из массива амплитуд:
3.1. Умножение массива амплитуд на оконную функцию;
3.2. Преобразование Фурье.
4. Фильтрация массива спектральных значений.
5. Расчёт массива громкости (модуль от спектральных значений).
6. Нормализация громкости:
6.1. Поиск максимальной громкости (maxVolume).
6.2. Нормализация (volume[i] = volume[i]/maxVolume);
7. Квалификация спектра как информативного или зашумлённого.
8. Инициализация 3 пиков громкости и соответствующих им частот.
9. Поиск фундаментальной частоты.
9.1. Если один из 2 главных пиков громкости ниже некоторого порога громкости, то выбрать частоту другого главного пика, выход из алгоритма.
9.2. Поиск фундаментальной частоты по 2 главным пикам согласно известным пропорциям 1/2, 1/3, 2/3, 2/5, 3/4, 3/5.
9.3. Коррекция фундаментальной частоты по 3-му пику.

Мне нужно пробежаться по некоторым пунктам, чтобы удостовериться в валидности алгоритма.
Пункты 1, 2, 6, 8 и 9 вопросов не вызывают.

Самые горячие вопросы сейчас такие:
1. Правильно ли считается громкость? Массив спектральных значений имеет следующую структуру: в ячейках spectre[2*i] лежат коэффициэнты при косинусах, в spectre[2*i+1] лежат коэффициэнты при синусах. Соответственно, громкость считаю так: volume[i] = sqrt(sqr(spectre[2*i])+sqr(spectre[2*i+1])); Может, таким методом не громкость считается, а что-то другое?
2. Как рассчитать шум? И как рассчитать порог шума (если его превысить, то выход из алгоритма)? Сейчас шум рассчитывается как среднее арифметическое по массиву громкости.

Сообщение отредактировал QuickNick - Jul 1 2011, 07:00

[SPACUM]

Мне нужно улучшить алгоритм поиска основной частоты

1.Во первыых терминология. До БПФ массив называется сигнал или выборка и алгоритм получения выборки не имеет отношения к заданным вопросам. После БПФ массив называется спектр, он комплексный Re и Im. После взятия корня из суммы квадратов получается массив амплитуд гармоник Фурье.
2.Термин "громкость" уже обсуждался и оказалось громкость = максимальная амплитуда.
3.На тему "поиск основной частоты" опубликовано множество работ и продолжает публиковаться. Это означает, что удовлетворяющий всех алгоритм не найден. Но существует множество приближений от самых простых приведенных в работах Куцукоса(Peter Kootsookos) до весьма сложных (Quinn, M. Macleod). Так что вполне можете придумать свой.
Можете поискать по словам "estimation" "estimation frequency" "funamental frequency" там много.
4.Рекомендую попробовать такой метод: Окно Гаусса->БПФ->Определить основные вершины->Уточнить вершины по трем точкам.
5.В измерении шума всегда используют критерий мощности. Возвести каждое значение сигнала в квадрат просуммировать разделить на количество и извлечь корень. Или возвести каждое значение амплитуды спектра в квадрат просуммировать и извлечь корень. (Делить не надо). Эти значения должны быть равны.

[QuickNick]

1.Во первыых терминология. До БПФ массив называется сигнал или выборка и алгоритм получения выборки не имеет отношения к заданным вопросам. После БПФ массив называется спектр, он комплексный Re и Im. После взятия корня из суммы квадратов получается массив амплитуд гармоник Фурье.
2.Термин "громкость" уже обсуждался и оказалось громкость = максимальная амплитуда.

Ага, спасибо, пока ещё "плаваю" в терминологии.

5.В измерении шума всегда используют критерий мощности. Возвести каждое значение сигнала в квадрат просуммировать разделить на количество и извлечь корень. Или возвести каждое значение амплитуды спектра в квадрат просуммировать и извлечь корень. (Делить не надо). Эти значения должны быть равны.

То есть, "должны"?
Они обязательно будут равны или же "если равны, то ..."?
Допустим, я получил величину "корень из суммы квадратов амплитуд спектра". У неё какая размерность?
Если её поделить на максимальную амплитуду, то это можно будет использовать как некоторый порог шума?

...бред какой-то )))... Вам нужно удостовериться, что главный пик это голос а потому имеет "ощутимую" 2-ую и третью гармоники...

В комментах предыдущих программист так написал:

Код

// If the volume is one of the 2-harmonics with maximum volume less than the threshold SIGNAL_LEVEL, it makes no sense to use to search for
// other harmonics of the fundamental frequency, because they are very weak, and perhaps are simply noise. Return as the fundamental frequency
// frequency with the maximum volume

[DRUID3]

В комментах предыдущих программист так написал:

Код

// If the volume is one of the 2-harmonics with maximum volume less than the threshold SIGNAL_LEVEL, it makes no sense to use to search for
// other harmonics of the fundamental frequency, because they are very weak, and perhaps are simply noise. Return as the fundamental frequency
// frequency with the maximum volume

...что переводится так - "если высшие гармоники ниже шумового порога нет смысла их искать, они "потоплены" в шуме"... Вы же сделали совсем странные выводы...

Кстати, SPACUM, дает Вам не очень грамотные советы(относительно шума) - он, видимо, программист от матлаба. Хотя может Вам тоже прога для матлаба нужна?

[QuickNick]

...что переводится так - "если высшие гармоники ниже шумового порога нет смысла их искать, они "потоплены" в шуме"... Вы же сделали совсем странные выводы.

Таков алгоритм предыдущего программиста. Повторюсь, всё, что я изменил в нём - это сделал настраиваемость оконных функций и фильтров (в предыдущей версии они были "железно" прописаны).

Кстати, SPACUM, дает Вам не очень грамотные советы(относительно шума) - он, видимо, программист от матлаба. Хотя может Вам тоже прога для матлаба нужна?

На Яве пишу.

[DRUID3]

Таков алгоритм предыдущего программиста.

Педыдущего если иди сверху? А если серьезнее - не знаю какой там алгоритм, но ваши пункты мало соответствуют тем
строчкам на английском языке. что Вы привели. Потому - не знаю что и думать.

Повторюсь, всё, что я изменил в нём - это сделал настраиваемость оконных функций и фильтров (в предыдущей версии они
были "железно" прописаны).

...а я бы оконные функции выбросил. Размен улучшения разрешения по частоте на уменьшение разрешения по времени(по-
сути потерю "энергетики") не оправдывает себя для поиска основного тона. И метрология здесь ни к чему, мало того с
окнами она ложная.

На Яве пишу.

Гы... Прикол. Не забудьте еще частоту дискретизации 8000 выставить, а то наваяете... А какая "ява"? "Наша" или
"андроидная"?

Я бы делал так:

1) скользящее среднее проверяет наличие звуковой активности(т.к. блок>>фонемы).
2) автокоррелятор проверяет шумовой ли сигнал (единственный пик на 0-е) или колебательный(периодическая
автокорреляционная функция). Ну и находим основную частоту - 1/период_автокорреляционной_функции .
))))))) . Но этот путь жрет много ресурсов. Можно автокорреляцию аппроксимировать конечно. Шагать непосэмплово. Но
это трудный путь )))...

...тогда так:
1) скользящее среднее проверяет наличие звуковой активности(т.к. блок>>фонемы).
2) если есть активность, то RealFFT начиная с момента активности.
3) Получение быстрой оценки физического спектра - Р.Лайонс ст.475
4) Получить массив указателей на массив значений физического спектра. Ну не знаю как там в Java придумаете. Нужно
чтобы и номера отсчетов первоначального спектра остались и отсортированный их массив был. Ну а далее техническое
творчество. Если наибольший спектральный компонент не на нулевой частоте и содержит ощутимую часть(вот этот бы порог
сделать изменяемым-подстраиваемым) от общей энергии (сумма всех компонент физического спектра) то осталось проверить
наличие 2-ой и третьей гармоники. Если такой компоненты нет - методом исключения - сигнал шумоподобный. Экономим кучу нервов и MIPS'ов .

2SPACUM:
Вы уж не серчайте, что я так о Вас - мол, программист от матлаба. Просто бывают и такие. Кому математический формализм важнее физической реализации. И 100% и такие нужны... Но... Как видите в предложенной мной схеме не нужно мерять шум возводя отсчеты в квадрат - экономим сотни MIPS'ов но не противоречим здравому смыслу.

[QuickNick]

Гы... Прикол. Не забудьте еще частоту дискретизации 8000 выставить, а то наваяете... А какая "ява"? "Наша" или
"андроидная"?

А на ПК какой библиотекой для записи звуков нужно пользоваться?
Хотелось бы приложение сделать, в котором рисовались бы промежуточные данные обработки сигнала. То есть, к компу подключается микрофон, и в реальном времени рисуются графики, спектры и т.д.

[SPACUM]

А на ПК какой библиотекой для записи звуков нужно пользоваться?
Хотелось бы приложение сделать, в котором рисовались бы промежуточные данные обработки сигнала. То есть, к компу подключается микрофон, и в реальном времени рисуются графики, спектры и т.д.

Пробовал только на билдере и давно, для него: http://edn.embarcadero.com/article/28332.
(!На сррbuilder.ru не заходите - это чистый вирус!, сворует все, можете врагу подсунуть)
Разные методы это чтобы уменьшить выборку и работать в реальном времени. Если между частотами меньше 5 бин все равно плохо получается.
А максимального правдоподобия - это так: Находите Вашим методом 3 частоты-амплитуды-фазы, всего 9 параметров и делаете оптимизацию по минимуму среднеквадратичного отклонения от исходной выборки. Результат превосходный, но слишком долго.

[QuickNick]

Пробовал только на билдере и давно, для него: http://edn.embarcadero.com/article/28332.
(!На сррbuilder.ru не заходите - это чистый вирус!, сворует все, можете врагу подсунуть)

Не, мне для Java нужно.
Но всё равно спасибо - когда-нибудь и C++ займусь.

Разные методы это чтобы уменьшить выборку и работать в реальном времени. Если между частотами меньше 5 бин все равно плохо получается.
А максимального правдоподобия - это так: Находите Вашим методом 3 частоты-амплитуды-фазы, всего 9 параметров и делаете оптимизацию по минимуму среднеквадратичного отклонения от исходной выборки. Результат превосходный, но слишком долго.

Бин - это что такое?
И какой тогда у меня метод получается?

Сообщение отредактировал QuickNick - Jul 6 2011, 06:38

[SPACUM]

Бин - это что такое?

Bin - расстояние между двумя ближайшими гармониками Фурье(1 / длина выборки).
А метод слишком грубый.
Например берем выборку длиной 1с. При этом 1бин = 1Гц. Подаем смесь трех кратных частот: 10.5бин - 2.0В, 21бин - 3В и 31.5бин .5В.
Получаем три максимальных гармоники 10бин - 1.41В, 11бин - 1.41В и 21бин - 3В, а 31.5бин в расчетах уже не участвует. И ничего кратного!

Сообщение отредактировал SPACUM - Jul 6 2011, 12:21

[QuickNick]

Bin - расстояние между двумя ближайшими гармониками Фурье(1 / длина выборки).
А метод слишком грубый.
Например берем выборку длиной 1с. При этом 1бин = 1Гц. Подаем смесь трех кратных частот: 10.5бин - 2.0В, 21бин - 3В и 31.5бин .5В.
Получаем три максимальных гармоники 10бин - 1.41В, 11бин - 1.41В и 21бин - 3В, а 31.5бин в расчетах уже не участвует. И ничего кратного!

А. Понял. Об этом я не писал, но после нахождения пиков мы уточняем частоту:

Код

    protected final double accuracyFreq(int index) {
        double res;
        if (volume[index + 1] > volume[index - 1]) {
            res = (index * volume[index] + (index + 1) * volume[index + 1]) / (volume[index + 1] + volume[index]);
        } else {
            res = (index * volume[index] + (index - 1) * volume[index - 1]) / (volume[index - 1] + volume[index]);
        }
        return res * sampleRate / signalLength;
    }

Насколько я помню университетскую программу, предыдущий программист здесь применил интерполирование.

Сообщение отредактировал QuickNick - Jul 6 2011, 14:00

[SPACUM]

А. Понял. Об этом я не писал, но после нахождения пиков мы уточняем частоту:

Код

    protected final double accuracyFreq(int index) {
        double res;
        if (volume[index + 1] > volume[index - 1]) {
            res = (index * volume[index] + (index + 1) * volume[index + 1]) / (volume[index + 1] + volume[index]);
        } else {
            res = (index * volume[index] + (index - 1) * volume[index - 1]) / (volume[index - 1] + volume[index]);
        }
        return res * sampleRate / signalLength;
    }

Насколько я помню университетскую программу, предыдущий программист здесь применил интерполирование.

Это не интерполяция, а одна из грубых формул Куцукоса.
Мой пример, должно получиться ~10.5бин.
res =(10 * 1.41 + 11 * 1.41) / (1.41 + 1.41) = 10.5бин
На самом деле если частота точно посередине между бинами значения справа и слева немного отличаются.
Если плавно менять частоту в этой точке будет ступенька, но для Вашего случая точность достаточная, более точные методы значительно сложнее.

Сообщение отредактировал SPACUM - Jul 6 2011, 16:23

[QuickNick]

res =(10 * 1.41 + 11 * 1.41) / (1.41 + 1.41) = 10.5бин
На самом деле если частота точно посередине между бинами значения справа и слева немного отличаются.
Если плавно менять частоту в этой точке будет ступенька, но для Вашего случая точность достаточная, более точные методы значительно сложнее.

Вот и нарвались на ступеньку
Так что теперь нужны более точные методы.

С одной стороны, ясно, что это будет задача интерполяции, но, с другой стороны, смущает неинъективность отображения в районе пика.
Какой математический инструмент на вооружение взять?

----Update----
Появилась следующая идея: построить интерполяционный многочлен по 3 точкам (freq[index-1], freq[index], freq[index+1]) и пробежаться с шагом 0.1 по нему. Найти максимум, запомнить значение.

Сообщение отредактировал QuickNick - Jul 14 2011, 07:26

[SPACUM]

Появилась следующая идея: построить интерполяционный многочлен по 3 точкам (freq[index-1], freq[index], freq[index+1]) и пробежаться с шагом 0.1 по нему. Найти максимум, запомнить значение.

1.Если 3 точки, то есть формула, еще в школе изучали. Если не помните - могу напомнить.
2.У меня по трем точкам получалось лучше с окном Гаусса и если частоты раздвинуты больше, чем на 5 бин.
3.Либо без окна и Куцукос, либо с окном и интерполяция.
4.Явление Гиббса никакого отношения к Вашим проблемам не имеет.
5.Пора уже написать проверочную программу и сравнивать точности. Примерно так:
Задать три случайные частоты->если какая-нибудь меньше 7бин или больше 80% рабочего диапазона или разность каких-нибудь частот менее 5бин то переделать-> добавить немного шума->определить 3 частоты любым Вашим методом и определить максимальную погрешность определения частот. Другим способом проверить метод невозможно.
6.Если ограничения не устраивают, то следует ввести расчет процента сбоев. Лучше сначала определить точность с этими ограничениями.

[QuickNick]

1.Если 3 точки, то есть формула, еще в школе изучали. Если не помните - могу напомнить.

Запамятовал
Но сейчас в другую крайность влетел.
Если Куцукос "ступеньку" выбрасывал, то сейчас подряд идут значения, отличающиеся всего на 1-2 десятых, но так как они сменяются каждые 200 мс, то прибор опять же остаётся неудобным.
Я сейчас сижу над задачкой увеличения стабильности. Период взятия данных не должен превышать 250 мс (лучше всего 200 мс), чтобы показывал точно и в то же время на мелочь не реагировал.

Опробую сейчас с различным окнами. Хотя не исключено, из малого опыта опять не в ту степь пойду.

P.S.: спасибо вам большое!

Сообщение отредактировал QuickNick - Jul 15 2011, 06:44