Я, как (к сожалению) далекий от теории человек, делал проверку метода экспериментально на модели: задавал чистую синусоиду и вычислял разницу между этой подаваемой частотой и той, что мне метод насчитывал. В цикле для разных частот, потом на график. Если зависимость немонотонная- то детально обследовал внутри этой немонотонности (у меня, например, достаточно между двумя бинами ДПФ исследовать). Далее то же, но с использованием именно тех типов данных и размерностей, которые у Вас реально применимы в железе, не все и не везде double имеют. Если на этом этапе не нравится- то надо что-то в консерватории менять.
Ну а далее- добавляете мусору в исходный сигнал по вкусу- хоть гармоники, хоть белый шум, хоть еще что-нибудь, присущее вашей физике исходного сигнала.

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

Прежде чем экспериментировать и бросаться в пучину расчётов, рекомендую рассчитать два массива данных для интересующей длительности наблюдения и сравнить их.
1. Допустим, речь о частоте Fи
2. Желаемая точность dFи
3. Частота дискретизации Fд
4. Разрядность АЦП - N
5. Период наблюдения - T.
Таким образом, имеем T*Fд выборок, квантованных (в лучшем случае) по 2^(N-1) уровням...
Вычисляем 2 идеальных массива для частоты Fи и Fи + dFи и сравниваем их. В ряде случаев эти массивы (при чрезмерно оптимистичном подходе ) будут просто идентичными, тогда и копья ломать не стоит.
Если идти чуть дальше, то можно соотнести невязку этих двух массивов и одним из них с величиной, например, сигнал/шум реального АЦП. Это тоже покажет вам примерные идеальные возможности вашей системы.

Ruslan1, подскажите пожалуйста адреса ресурсов, где вы нашли функции по вычислению окна Гаусса и гаусовской интерполяции.
К примеру, на сайте http://www.dsplib.ru/content/winadd/win.html указана немного другая формула по вычислению коэффициентов окна Гаусса.

Просто я попробовал применить ваш алгоритм на С++ на реальном сигнале от гитары, так он выдает всегда результат интерполяции на 0,5 больше чем исходная частота.
В математике я слаб, но ваш метод интерполяции мне показался уж очень простым

http://www.nicholson.com/rhn/dsp.html - там есть подборка ресурсов на тему уточнения FFT
http://www.mathworks.com/help/toolbox/signal/gausswin.html


Да каких только формул я не видел, пока искал
Я уже писал, что не уверен в том, что мой вариант верный. Но матлаб для меня это авторитет, а у него внутренняя функция так же работает как я считаю.


0.5 Гц? странно, о такой разнице я еще не слышал. Можно исходники и окна и интерполяции в студию? И файл с исходным оцифрованным сигналом, очень интересно.....
Если на 0.05 Гц- это уже ближе к телу....

Сначала я в одном источнике прочитал что точнее будет, если логарифмическую кривую аппроксимировать параболой, а потом совсем в другом месте прочитал что это и есть гауссовская интерполяция.

Правильное Гауссовское окно считается так:
sqrt(sqrt(PI)) / sqrt(sigma) * sqrt(FSIZE) * exp (-(x-FSIZE/2.0)*(x-FSIZE/2.0)/2.0/sigma/sigma);
PI - это число пи, FSIZE - размер Фурье в точках, sigma - параметр, задающий ширину окна.
При таком нормировочном коэффициенте накладывание окна не приводит к дополнительному коэффициенту в амплитудах.
По поводу гауссовской интерполяции - все правильно, но полезно вводить весовые коэффициенты, пропорциональные квадрату амплитуды сигнала - улучшает точность, т.к. шум при малом сигнале меньше влияет на определение положения вершины.

Извиняюсь за очередной глупый вопрос: а как выбрать ширину окна (sigma)?
Например, хочу задавить третью гармонику, которая в худшем случае находится на расстоянии 105 бинов, если рассматривать частотную область. Это значит, что мне нужно взять окно, которое в временной области задавит находящийся на расстоянии 105 семплов от ценрального отчета (FSIZE/2) ? Или все совсем не так?
У меня получилось что при сигма=43 имею подавление 60 dB во временной области.

Или достаточно, чтобы к границам окна наблюдения (отчеты 0 и FSIZE) упало до нуля ( ну или -60 dB хотя бы) ?

Ширина выбирается в зависимости от задачи. Мне нужно было, чтобы вклад второй гармоники в амплитуду основного сигнала не превышал 1е-8, ну я и подбирал сигму так, чтобы посередине между основным тоном и второй гармоникой был провал не хуже 1е-8. Для определения только частот можно, думаю, ограничиться разделением на уровне 1е-6.
По поводу соотношения окна в частотной и временной областях - что-то там похожее должно быть, но я не исследовал точно. Может какой-то коэффициент вылезти. Лучше смотреть сразу в спектре.

Ruslan1, извиняюсь, при копи-пасте потерял еденичку в формуле интерполяции. Пока дебаггером не прошел, не замечал.
За ссылки - большое спасибо, немного подразобрался с теорией.

да пожалуйста.



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

Шкалу вертикальную сделайте линейную, а не логарифмическую - тоже ничего не увидите. У них же не децибелы, а амплитуда.

Строго говоря, приводит, т.к. коэффициент у вас рассчитывается для всей прямой, а не для ограниченного носителя.

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

Вопросы о разрешающей способности измерения частоты за доли секунды с разрешением меньше Герца можно посмотреть здесь, очень невероятные вещи узнаете!
http://electronix.ru/forum/index.php?showt...mp;#entry703059

Для струны (струнные датчики растяжения) есть еще две проблемы, мы на них наезжали.

1. Частота колебаний зависит от их амплитуды. Грубо говоря, если струна колеблется с большей амплитудой - она больше растянута и ее жесткость выше. Единственный выход - минимизировать амплитуду ((или стабилизировать или выбирать участок с определенным ее значением).

2. В струне, особенно не очень корректно закрепленной, присутствуют колебания на нескольких очень близких частотах. Это связано, например с тем, что струна имеет чуть разные длины и жесткости для вертикального и горизонтального колебания, например на гитаре. А горизонтальная и вертикальная моды переходят одна в другую, в результате появляется медленная модуляция амплитуды на фоне общего затухания сигнала.

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

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

Интересно, спасибо. А еще, мне всегда казалось, что струна выдает не чистую синусоиду, а все-таки с гармониками, как минимум с четными, то есть у струны отдельно есть колебания ее половины, четверти...