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

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


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

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

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 - Jul 6 2011, 14:00

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

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

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 по нему. Найти максимум, запомнить значение.

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

Явление гиббса еще никто не отменял. Методы борьбы с ним широко известны. Возьмите окно кайзера с параметром ~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.: спасибо вам большое!

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