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

Мне нужно улучшить алгоритм поиска основной частоты по методу максимального правдоподобия. Был пример у меня, как делать, я его модифицировал в нескольких местах (добавил настраиваемость оконной функции и фильтра). Возможно, некоторые места вам покажутся совершенно абсурдными - но я пока новичок в 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. Как рассчитать шум? И как рассчитать порог шума (если его превысить, то выход из алгоритма)? Сейчас шум рассчитывается как среднее арифметическое по массиву громкости.

Описано много всего. Зачем все это делается - не понятно. Считаем, что требуется все-таки выполнить оценку частоты и амплитуды гармонического сигнала. Требования к качеству оценки не предъявляются. Таким образом, решение - накопить достаточное количество отсчетов сигнала и взять БПФ от них. Получите оценку частоты, начальной фазы и амплитуды.

Нужно в реальном времени давать частоту сигнала.
То, что БПФ надо брать - это понятно.
У меня вопросы только по конкретным шагам - насколько они валидны.

Собственно, я нашёл решение по шагу 7. Посчитать среднее арифметическое по громкости в полосе фильтрации, сравнить его с пиком громкости. Если шум составляет >10% от всей информации, то далее последовательность не обрабатывать.
Для этого я шаги 7 и 8 местами поменял.

Алгоритм работает, всем спасибо за сочувствие!

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

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


...в измерении - возможно, но не в оценке - уж точно... Именно математический формализм отъедает основные "мипсы" в прожектах новичков ...

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


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


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

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

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

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


На Яве пишу.

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


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


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

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

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

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

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

А чего серчать? Каждый может ошибаться. Матлабом не пользуюсь совсем. Все алгоритмы о которых пишу проверил только на микроконтроллерах и на реальных сигналах. Люблю 32-битную целочисленную математику. Всякая отсебятина требует тщательной проверки на синтетических и реальных сигналах. В чудеса я верю, но это должны быть надежно проверенные чудеса.
Описанный мной метод показался мне надежным и устойчивым, только 50Гц и 100Гц без сигнала обязательно влезут. А качество сигнала легко вычислить зная суммарную энергию найденных трех гармоник и суммарную энергию всего сигнала.

Еще раз прошу прощения, что так Вас обозвал ...


...дело не в том, что Вы ошиблись - нет это была не ошибка, а в том что немного не туда увели беседу... И очень было похоже на "акадЭмика" который вобщем-то в теме, но в виду специфики своего опыта дает непрактичные советы... Понимаю - я был неправ.


Тем более, думаю, Вы оцените, что своими, возможно иногда и бестактными советами я таки да и облегчаю жизнь разработчикам ...


...так все оно отсебятина... Ряды Фурье так невзлюбили поначалу - "отсебятина какая-то" - и кто ученные в субстанцию моченные мирового уровня, а сейчас рядового инженеришку за ухи не оттянуть...


Так выкинуть их из FFT... обнулить...

До меня программист писал - как он объяснил, таким методом он в жизни струны настраивает. Может, я его неправильно понял, не знаю - просто музыкального слуха у меня нет, гитару настроить не смогу никогда.
Я потестил - впечатление, что действительно не нужен пункт с проверкой вершин относительно "сигнального уровня".


Тестеры разберутся, что оптимальнее. Нужно будет - поставим прямоугольное окно.


Андроидная Ява, да.
С 8000 уже столкнулся - решил проверить, что будет, если увеличить частоту дискретизации вдвое.
А действительно, почему так?
Что произошло?


Хух Столько терминов, что у ученика перегрузка мозга произошла

Спасибо!

У Вас английский "внутрикорпоратив"?


...им не пофиг?


Ха... а что там есть из звукового API?


Х.з. ... а если серьезнее - я вообще о том, что нужно выставлять 8000 во всех приложениях связанных с голосом иначе оверхед - ну там 44100 например. А не заработало у Вас скорее всего из-за того, что звуковушка там 16000 не поддерживает.


Скользящее среднее тотже Лайонс ст.398. RealFFT кажется есть в API "андроида"... "квиксорт" тоже должен быть. А о автокорреляторе просто пока забудьте.


Стоп, а деньги где?

P.S.: ...простой проект по захвату звука в "андроиде"... Вдруг кому пригодиЦЦо...

1. Комменты он писал на английском - не знаю зачем.
2. Они скажут, при каких параметрах тюнер лучше работает - их уже намертво пропишем. Чтобы легко всё настраивалось, я в виде фабричных методов всё организовал.
3. Не в курсе.
4. Не заработало из-за приведения типов - я уже разобрался.

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

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

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


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

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

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

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

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

1)Есть у вас пример на руках?
2)Хотелось бы сделать в текущем алгоритме оптимизацию по фильтрации (подавать на вход предполагаемый диапазон частот). Но из нашего алгоритма (как он называется, кстати? Метод обертонов?) получается, что если мы в фильтр подадим правую границу диапазона, то будут сильно глушится 2-я, 3-я и далее гармоники. Поэтому сейчас приходится идти на ухищрения: подавать утроенную правую границу диапазона. Есть ли алгоритм, который бы позволил избежать таких ухищрений? Конкретно, можно ли алгоритм максимального правдоподобия без таких вывертов писать?

1. Я программирую для архитектуры АРМ, сейчас для Кортекса. Если есть исходник, то он либо применен в оригинальном приборе, либо это задел.
В обоих случаях раздавать глупо.
Два алгоритма описал достаточно, можно сразу программировать.
В случае оптимизации: амплитуды надо нормировать к найденным значениям, смещение фазы от найденной к ПИ, смещение частоты от найденной к 1 бин. Можно применить покоординатный спуск или градиентный по каждой триаде. Пробовал. Очень долго. Не рекомендую. Просто обьяснил смысл термина.
2. А если нужен железный пример на руках, приезжайте на следующую конференцию NXP в Москве, на две прибор притаскивал и на следующую притащу. Там есть поиск четырех реальных частот в шумах, разумеется не методом оптимизации. Интерес низкий.
3. Если Вы применяете метод Куцукоса, никаких окон делать нельзя, с окнами метод показывает значения точно посередине между ближайшими двумя гармониками Фурье и ничего не уточняет. А вообще выше 7-й гармоники Фурье и при малых шумах Куцукос весьма точен.
4. Не понимаю зачем фильтрация? Спектр получен и Вы можете в нем выбрать интересующий Вас участок.
Вот аналоговый фильтр перед АЦП надежно не пропускающий частоты выше или равные половине частоты дискретизации обязательно должен быть.
Например если частота дискретизации 8кГц, то фильтр должен быть с частотой среза 3.5 кГц чтобы частоты 4кГц и выше были полностью зарезаны, иначе все недорезанные частоты выше 4кГц будут присутствовать в диапазоне рабочих частот и никакой программой их оттуда не убрать.
5. Поиском "fundamental frequency" не занимаюсь, речь шла только об определении реальных гармоник вместо гармоник Фурье.

2 SPACUM: спасибо большое за помощь!

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

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

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

Явление гиббса еще никто не отменял. Методы борьбы с ним широко известны. Возьмите окно кайзера с параметром ~30. Будет Вам равенство значений соседних точек, если частота попадает точно в середину между бинами.

Уже расправился - пока использую интерполяцию Ньютона по 5 точкам, затем с шагом по 0.1 Гц прохожу по интервалу (он получается длиной примерно 3.9 Гц), нахожу максимум.
Возможно, есть какой-то хитрый способ интерполяции наоборот, когда отображение не инъективно, но за полдня найти его не успел.

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

Это и есть явление гиббса.

Небольшой разрыв в графике определенной частоты от исходной при использовании формулы Куцукоса. Если применить любое окно, то погрешность равна 1/2бин вместо долей процента для 100-й гармоники. Куцукос работает только с прямоугольным окном как и было в исходной программе. Рекомендация не интересна.
Также не вижу нигде 18% предсказанных Гиббсом, а их всегда видно.

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

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

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

Если y0 - значение посередине для частоты x0, y1 - значение слева для x0 - dx, y2 - значение справа для x0 + dx, dx - шаг по частоте, то
x = x0 + dx * (y2 - y1)/(2 * y0 - y1 - y2) / 2;
Это экстремум параболы по трем точкам.

Точно, спасибо!

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

Собственно, так и делаем. Плюс ещё сравниваем поведение с другими детекторами.

А стабильности я добился с помощью очереди, которая в TuxGuitar реализована (org.herac.tuxguitar.gui.tools.custom.tuner.TGTunerQueue). Вернее, её взял за основу.