Что-то я не знаю приборов измеряющих арктангенс произвольного аналогового сигнала, может просветите?

Кроме того, арктангенс - разрывная функция и брать производную от него...

Несерьезно как-то, лучше книжку почитайте.

Когда-то в 1993г. такую задачу решал для измерения токов в рельсах метро на контроллере 8051 на частотах 50, 75, 125 ... Гц.
Решалась просто: сигнал умножался на косинус, синус 50 Гц и накапливался как скользящее среднее. Результат - как фильтрация узким фильтром. Если сигнал строго =50 Гц. то выходной сигнал - постоянная, если 49 Гц - вых. сигнал 1 Гц, если 51 Гц - сигнал 1 Гц, но крутящийся в противоположную сторону. И т.п. Соответственно, если на выходе синусоида с периодом 100 сек. то сигнал 50,01 Гц.

дак это надо тогда 100 секунд ждать, пока хоть 1 период наловишь? А чтобы принять надёжное решение, нужен ещё и не один период? Видимо хотелось бы побыстрее...

Просвещайтесь на здоровье

А чего это они в формулах (4) и (7) после второго знака равенства "arctan" выкинули?

Арктангенс равен аргументу комплексного числа, записанного в (3) как две отдельные квадратуры, обычное приравнивание. Если кому-то интересно, то вот о процессе угловой модуляции. Тогда всё сразу становится понятнее )

Да читали мы все это. И даже экзамены сдавали по радиолокации.. Просто очень давно..

И где же тут прямое измерение "мгновенной" частоты? После фильтров ФНЧ?
Вы просто подтвердили мои слова:

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

И, как обычно, вы перепутали реально существующую и оцениваемую по достаточно длинной выборке разность фаз двух сигналов с мифической фазой произвольного ОДИНОЧНОГО аналогового сигнала.

Измеряется не фаза, а разность фаз, и за несколько периодов.

И, если вы прочитали всю приведенную ссылку, в ней, опять же, подтверждают мои слова о разрывности дифференцируемого вами арктангенса

Поэтому, в реальной обработке, арктангенс старательно исключают.

Расскажите мне, как определить фазу сигнала по одному отсчету АЦП? Ну ладно, учитывая нарисованные вами производные - пусть будет целых три последовательных отсчета с интервалом 1 мксек: -1В; 0В; +1В.
Ну, хотя бы, любимую "мгновенную" частоту определите?

Не получается? Тогда прекратите увиливать, лучше читайте книжки и думайте.

Нет ничего проще: На основании именно этих полученных от АЦП отсчетов можно смело сказать, что на его входе присутствует синусоидальный сигнал с частотой 250 kHz, c амплитудой 1В. Центральному отсчету соответствует фаза 0. Мы считаем, что амплитуда сигнала нормирована под полный "раскрыв" АЦП +/-1В.

Прошу прощения за нескромный вопрос. А в каком месте atan(x) разрывается?

И прочитал и в своё время сделал в железке, не поверите
Его там исключают потому что далее делается дифференцирование и можно вместо d/dt(arg(y(t))) рассчитывать конечное выражение и собрать эквивалентую схему. Никто не запрещает сделать схему "в лоб" с использованием арктангенса + unwrap и всё будет работать (проверено на практике, работает, alas!).
Да и где ж разрыв?
Нажмите для просмотра прикрепленного файла


Приехали. Преобразование координат и демодуляцию PSK-N созвездий наверное тоже без арктангеса делают, ага.


А мне не нужно работать на уровне АЦП. И ТС'у никто не предлагал работать на уровне АЦП и по одному отчёту. Есть сигнал после гетеродина и ФНЧ, и он комплексный. Для него понятие фазы определено строго математически и применимо к каждому символу. А раз определено понятие фазы, то значит можно оценить мгновенную частоту. Точность такой оценки будет зависеть от сигнал-шума. Процесс демодуляции ЧМ сигналов привёл как простой и распространённый пример.

Что ещё Вы хотите тут опровергать?

Наверное, подразумевалась неоднозначность.

Можно вполне строго доказать невозможность измерения "мгновенной частоты" на практике.

Для этого нужно воспользоваться формулой для вычисления частоты по трем значениям гармонической функции S(t) из соседней темы:

f = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

Это абсолютно точная формула, позволяющая получить точное значение частоты "f" при условии,

что нам известны точные значения гармонической функции S(t) в моменты времени t₁,t₂ и t₃:

S₁ == S(t₁),
S₂ == S(t₂),
S₃ == S(t₃),

где:

Δt == t₃ - t₂ = t₂ - t₁,

и при этом:

f*(t₃ - t₁) < 1.

Чтобы строго доказать невозможность измерения "мгновенной частоты" на практике, необходимо вычислить полный дифференциал функции трех переменных:

f(S₁,S₂,S₃) = [1/[pi*(t₃ - t₁)]]*arccos[(S₁ + S₃)/(2*S₂)].

А именно:

df == (δf/δS₁)*dS₁ + (δf/δS₂)*dS₂ + (δf/δS₃)*dS₃,

где частные производные функции "f" равны:

δf/δS₁ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

δf/δS₂ = +[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(S₁+S₃)/(2*S₂²)].

δf/δS₃ = -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[1/(2*S₂)],

Соответственно, полный дифференциал равен:

df == -[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*[(dS₁+dS₃-(S₁+S₃)*dS₂/S₂)/(2*S₂)].

Теперь мы можем найти погрешность измерения частоты "σ_F" связанную с погрешностью измерения сигнала "σ_S":

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[1-[(S₁ + S₃)/(2*S₂)]²]]*sqrt[2+(S₁+S₃)²/S₂²]/[2*abs(S₂)],

или:

σ_F = σ_S*[1/[pi*(t₃ - t₁)]]*[1/sqrt[(2*S₂)²-(S₁ + S₃)²]]*sqrt[2+(S₁+S₃)²/S₂²].

Знаменатель в этой формуле примечателен тем, что при стремлении t₃ к t₁ он стремится к нулю как const*(t₃-t₁)²

Как следствие, при стремлении t₃ к t₁ для любой погрешности "σ_S" измерения амплитуды сигнала S(t), погрешность измерения "мгновенной частоты" "σ_F" устремляется к бесконечности как σ_S/(t₃-t₁)²..

Смелый вы человек. Но с чего вы взяли, что это синус? Боженька на ушко шепнул? А может это линейно нарастающее напряжение? Или пила, трапеция и т.д. с произвольным периодом?
Или вообще непериодический сигнал?

Но, даже если предположить чистый синус, то с чего вы решили, что амплитуда сигнала равна 1В?
Вот набор пар амплитуда/частота который дает именно эти отсчеты. Ваш вариант - это только один из них.
Амплитуда ----------- Частота --------- Напряжение при t=1 мксек
1.000000E+00 -- 2.500000E+05 1 -- Ваш вариант
1.500000E+00 -- 1.161398E+05 1
2.250000E+00 -- 7.329944E+04 1
3.375000E+00 -- 4.787579E+04 1
5.062500E+00 -- 3.164613E+04 1
7.593750E+00 -- 2.101973E+04 1
1.139063E+01 -- 1.399046E+04 1
1.708594E+01 -- 9.320293E+03 1
2.562891E+01 -- 6.211555E+03 1
3.844336E+01 -- 4.140452E+03 1
5.766504E+01 -- 2.760129E+03 1
8.649756E+01 -- 1.840034E+03 1
и далее ...

А если взять частоты более 250 кГц (не совсем корректно, но чего только в жизни не бывает) - то снова получим неограниченный набор пар частота/амплитуда, котрые дадут те же отсчеты.

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


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


Действительно не совсем корректно выразился. Разрывность появляется при попытке определить именно фазовый угол (что, при грамотной обработке, практически не требуется) по соотношению Re и Im компонент применяя функцию арктангенса. Нужна, как минимум, atan2, но и с ней получаемый фазовый угол будет не линейно нарастающим (если есть постоянная разность частот), а скакать при переходе через +/-пи.


Это как минимум.

А тыкните носом в раздел анализа и теории оценок, где написано, что надо так считать. Просто интересно
В моём понимании f = [1/[pi*(t3 - t1)]]*arccos[(S1 + S3)/(2*S2)] при t3-t1->0 примет вид acos(1)/0, что даёт ноль в числителе и ноль в знаменателе, нормальное поведение для производной.
Кстати ваш вывод будет справедлив только для гармонического сигнала, но всё равно любопытно.

С другой стороны для идеального (нет шума) действительного гармонического сигнала можно написать такой простой скрипт:

который вычисляет мгновенную частоту по двум точкам (по одной точке можно оценить фазу, но для действительного сигнала она не имеет смысла, т.к. в сущности может быть любой в зависимости от выбора начальной точки; для частоты достаточно двух точек). Легко проверяется в матлабе.
Увеличивая частоту дискретизации, мы сокращаем время между двумя точками, с другой стороны разница напряжений в этих точках тоже уменьшается, т.о. мы приближаемся к определению мгновенной частоты. Если шума нет, то мгновенная частота будет определяться с точностью до fs, если добавить шум, то погрешность, которую он внесёт, превысит погрешность от fs (при fs >> искомой частоты). Погрешность такого эстиматора уже будет оцениваться с помощью CRB.


Увиливаете тут только вы. Если вы хотите, чтобы ваши агрументы были приняты, опишите это математически или, если удобнее скриптом. Мы не на уроке литературы. А за одним посмотрите, чем отличается неоднозначность обратных тригонометрических функций от понятия разрыва и непрерывности. Это довольно познавательно.

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

Если следовать этой вашей логике, то любая метрология, в основе которой лежит приницип причинности, является баловством и никчемной ахинеей: собрали установку, измерили частоту в розетке, записали результат, но как раз в этот момент напряжение могут отключить. Зачем суетились, и что делать с записями - не понятно. Это первое.

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

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

А существование понятия "мгновенная частота".. Угловое положение вращающегося диска можно измерить. Время тоже.

Читайте про преобразования случайной величины, Якобиан и проч.

Или: Г.Корн Т.Корн, стр. 561, Гл. 18.5-4, "Функции от случайных величин. Замена переменных."

Ну и имейте ввиду, что случайная величина "f" есть однозначная функция трех независимых случайных величин: f = f(S1,S2,S3)..

Для неопределенной амплитуды скрипт, естественно мертвый. 3 точки нужно минимум.

Да, точно. Ну и на практике никто так не делает (x1 = sqrt(1 - x.^2)), если нужно получить аналитический сигнал, то ставят Гильберта. Там никакой зависимости от амплитуды + класс входных сигналов не ограничен только гармоническими. Ну или гетеродин + ФНЧ - если диапазон измерения априори задан.

Единственная формула, которую сразу можно посчитать в C++ (IDE NetBeans, компилятор MinGW ). Два результата: для 16-ти разрядного АЦП и для идеального АЦП (см. результат в конце C++ кода). Для хорошего АЦП метод точный, т.е. если фильтрануть основную гармонику, то будет хорошо.

Другие способы (MUSIC, MLE - Метод максимального правдоподобия...), по крайней мере, на первый взгляд очень сложные. Какие-то для матлаба, какие-то в общем виде многоэтажные формулы. Кто такие формулы пишет? Шифротекст для потомков. Неужели один я такой тупой, что не понимаю.

"Сложные проблемы всегда имеют простые, легкие для понимания неправильные решения."

Вы никак не улучшите ваш сигнал оцифрованный 16-битным АЦП. При таком соотношении символьной и искомой частот (3200/50) аргумент функции acos -> 1, а значит сама функция - к 0, т.е. некоторому малому числу. Далее вы делите на 2pi*dt, другое малое число. С вычислительной точки зрения мало устойчивая задача. Любое незначительное искажение сигнала внесёт значительное искажение в оценку. Пропустив сигнал через АЦП вы добавили шум квантования, который внёс ошибку в ваш эстиматор. А теперь представте - в реальности при оцифровке стараются сделать так, чтобы оцифрованный аддитивный шум (тепловой шум, шум канала, кароче говоря неустронимый физический шум) был много больше шума квантования - условие необходимое для сохранения SNR на том уровне, который был до оцифровки. И посчитайте, какую ошибку внесёт он. Будете неприятно удивлены, т.к. даже незначительный шум даст полностью бесмыссленные значения частот. И даже узкополосный полосовой фильтр на 50 Гц вряд ли существенно улучшит результат.
Так что вам придётся применять какой либо из "сложных" методов.

Так Вы, оказывается, не поняли основной идеи предложенной формулы..

Я же там где-то раньше указывал, что для получения точного значения частоты должно выполняться условие: abs(S₁ + S₃) < abs(S₂).

Это требование фактически означает, что для вычисления "f" годятся не любые моменты времени t₁,t₂ и t₃, а только те,

в которых значения S₁ и S₃ находятся "вблизи" нуля и при этом значение S₂ "автоматически" оказывается вблизи максимума функции S(t).

Это примерно соответствует вычислению частоты "вручную", то есть когда мы просто измеряем моменты времени в которые гармонический сигнал S(t) пересекает ось абсцисс,

потом находим соответствующий этим моментам времени период гармоники "T" и уже затем вычисляем её частоту f = 1/T. Просто в предложенном способе не нужно точно вычислять моменты времени, когда S(t) = 0.

PS. А про шумы полностью согласен. Формула с arccos применима только для сферических функций в вакууме..

Формула сама интересная и, если сигнал будет чистый, то должна дать хороший результат.
Попробую сигнал с выхода АЦП обработать НЧ фильтром с КИХ и полосой 75Гц -3 дБ, 100 Гц -40 дБ, 64 коэффициента float.
А там можно разные способы попробовать. Что-то похожее советовал мне petrov.

Только вот не пойму, откуда бины возьмутся после фильтрации сигнала фильтром с частотой среза 75 Гц, вроде там должен быть один бин, т.е. 50 Гц. Или я что-то не понимаю. Или он имеет ввиду посчитать 1024 точечное ДПФ от сигнала на интервале 0,02 секунды и по тем бинам посчитать?

Эту задачу на форуме решали уже десятки раз..

См. например:

И там ещё в окрестностях есть что почитать..

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

Pridnya, если с частотой семплирования и разрешающей способностью всё понятно, то на каком минимально допустимом интервале времени (в периодах основной гармоники) Вы бы хотели измерить частоту?

Хотелось бы измерять за 1, максимум за 5 периодов. Не знаю пока сколько достаточно и на сколько достоверно будет при большом интервале, а вдуруг там какое качание частоты, поэтому хотелось бы за один период.

Делов то. Фурье на 8к. Самая яркая палка это она))). Если предварительно сузить полосу прикинув частоту и лисперсию с помощью какого-нибудь капчюр модуля МК. Можно сделать перенос на 0 вокруг, то можно понизить частоту дискретизации и Фурье станет менее требовательным к числу точек и скорости вычислений.

Почитайте внимательно тему.
Фурье 8192 вы собрались на мк сделать?
На частоте дискретизации, обозначенной ТС'ом в 3200 Гц, 8192 точки это 2560 мс. 1 - 5 периодов 50 Гц сигнала это 20 - 100 мс.

Накидайте ка на Си, а я код прогоню на ARM и оценим:
а) реализуемость;
б) время выполнения кода в целевом устройстве;
в) глубину бездны знаний.

Сколько можно уже фигней страдать. Вот сюда посмотрите база сигнала
А потом положите B=1, дельта F=0.01 и вычислите дельта t. Затем примите как данность, что Вы проектируете показометр и рассчитывайте значение частоты хоть по 3 точкам на период хоть до 10 знаков после запятой.

Δ t = B/Δf = 1/0.01 = 100 (секунд?). Не мойму, что-ли 100 секунд предлагается ждать? Подобный метод уже предлагали (если на выходе биение 0,01 Гц, то частота 49,99 Гц).

Сигнал: синусоида f=50 Гц со случайной начальной фазой. fd=44100 Гц.
М- число периодов синусоиды.

Не понятно, как вы получили результат (матлаб посчитал с помощью dll-ки или вы на С код написали)? Понимаю, что частота дискретизации для оцифровки звука. Хотелось бы попробовать ваш метод расчета на С или С++ при отклонении частоты.

Интервал 5 периодов сети 50 Гц - это 0,1 секунды. При частоте дискретизации 44100 Гц, получаем 4410 точек выборки на том интервале. А дальше что делать?

Это скорее для иллюстрации, что на 5 периодах погрешность получится 0.25Гц, а вам надо 0.01 :-)
Алгоритм сложный, в МК точно не поместится.
Советую не заморачиваться, поставить 1-бит АЦП (0 - компаратор) и счётчик.
На этом принципе все протонные магнитометры работают. (Погрешность 0.005 Гц, частота 2-3 кГц, время 1-2 сек)

Да мне хотя бы 0,25 Гц программно получить, а там я бы посмотрел что дальше делать. Все говорят "возьмите...", а вот сколько точек на каком интервале взять и как дальше, вот что главное.
Хотя бы с пяти периодов начнем (сигнал на интервале 0,1 секунды). Частота выборки пока 3200 гц. Сколько точек БПФ брать на этом интервале?
И правильно ли я понимаю, что я получу разрешение по частоте 10 Гц? 64точки*5периодов = 320 точек. 3200Гц/320 = 10 Гц. Так?


У меня использовался компаратор, модуль захват-сравнение, можно вытащить 0,01 Гц, с той же звуковой карты компа подаю и вижу 49.00 Гц, 50,00 Гц...Но мне интересно, как это сделать программно (без матлабов и прочих коммерческих продуктов для настольконго компа) для микроконтроллера.
Заодно и пойму, как получить высокое разрешение по частоте.

Да так. При данной частоте дискретизации и N=320 получится именно такое разрешение по частоте.
Весьма наивно полагать, что проведя параболу через 3 точки, близкие к соответствующему максиммуму спектра, можно получить хорошее разрешение по частоте.

Уже лучше, можно смоделировать: задать чистую синусоиду 50Гц на интервале 0,1 сек, а потом от этого интервала взять 1024-х точечное ДПФ для нескольких частот, близких к искомой, а там посмотреть какая кривая получится.

А вот как считать? В смысле спектр будет симметричный относительно искомой частоты (чтобы парабола получилась), т.е. отрицательные и положительные частоты? Или от 0 до частоты выборки Fs.

Библиотека CMSIS_DSP для ARM выдает результат от 0 до Fs, т.е. спектр симметричный относительно Fs/2.

Можно сразу задать частоты 49,50,51 Гц и найти значения спектра на этих частотах. Потом провести параболу и найти её максимум. Это будет первое приближение по частоте. Затем взять +-0.5 Гц от найденой частоты, снова найти значения спектра и снова проводить параболу, находить ее максимум, делить интервал на 2 (т.е. до 0.25 Гц) и т.д. Такая вот итерационная процедура... Но это долго получится.

Буду пробовать. Спасибо! Открытие для себя сделал: оказывается можно больше одного периода сети анализировать, мне доступно 16, 32 ,64, 128, 256, 512, 1024, 2048 и 4096-точечное БПФ (спектр посмотрю). А я почему-то думал, что БПФ нужно на интервале периода брать, т.е. на интервале 0,02 секунды 64, 128, 256 и т.д. точек.

PS: Смущает только, как оценить скорость изменения частоты при больших интервалах. Но это вопрос следующий.

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

До меня только сегодня дошло, как увеличить разрешение по частоте, используя преобразование Фурье и длинные интервалы, но это при условии, что частота основной гармоники не меняется. А вот что делать если она медленно меняется и нужно еще знать скорость изменения частоты? Вроде, интервал нужно уменьшать, но тогда разрешение по частоте уменьшится. Фурье не подходит, вроде. Может вейвлеты?

А суть идеи я понял (я не понимал, что вы считаете бинами, я их гармониками называл и думал, откуда там бины рядом с частотой основной гармоники, а там абстракция), и статья Эрика Якобсена у меня отложена, спасибо, вы мне её и давали.

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

http://www.ibm.com/developerworks/ru/libra...rl_1/index.html

А про "давно бы уже...", так есть ещё и другие вопросы.

Должно дойти, что нужна точность, а не разрешение.



А Ньютону и листка бумаги хватит. Что не понятно в алгоритме, комплексный тон не можете сгенерировать, не понятно как 1 бин ДПФ вычисляется, слишком сложна самая простая из возможных формула Якобсена?

Элементарно:
берем отсчеты АЦП и у каждого отсчета мнимую часть приравниваем нулю, далее обрабатываем комплексный входной сигнал библиотечными функциями CMSIS DSP.

Можно и без комплексного руками посчитать в тетрадке (для каждого бина ДПФ). Раньше я брал 64-х точечное ДПФ на интервале 0,02 сек и получал частоту 1 гармоники 50 Гц, второй 100 Гц и т.д. И не понял про какие бины возле 50-ти Гц вы мне говорили (откуда думаю там бины, если сигнал всего один, 50 Гц, если задать только один сигнал) .Теперь понятно.

На мк всю жизнь использовали рекурсивный алгоритм Герцеля, если нужно было в real-time считать отдельные частотные бины. Есть его модификация на тему скользящего Фурье преобразования (sdft), про это есть и на форуме и в гугле. Хотите отслеживать флуктуации частоты в окресностях 50 Гц - заменяете вычисление отдельных бинов dft на sdft, а далее - как вам советовали. Получаете скользящую модификацию алгоритма. Разрешение при этом будет зависеть от задержки преобразования (задержка sdft аналог количества точек dft), а точность - от сигнал-шума, нижнее значение которого вы нам кстати так и не озвучили (с этого нужно начинать проектирование любого эстиматора).

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


dft - ДПФ?
sdft - буква s что означает? Может sign (знак?) Google не знает.
https://www.google.ru/search?ie=UTF-8&h...&gws_rd=ssl

А отношение сигнал/шум большое, поэтому не критично.

dft=discrete fourier transform
s=sliding.
Вы заявили что хотите точность 0,01 Гц. Для такой точности нет некритичного SNR. Это большая точность. Тем более, если вы хотите отслеживать изменение частоты во времени -> статичный спектральный анализ уже отметается.


Вот тоже самое. dft и sdft имеют одинаковый результат каждые N отчётов. sdft позволяет отслеживать изменение комплексных значений выбранных частотных бинов во времени с задержкой N. Обрабатывая их соответствующим образом (интерполяция например), вы получите слежение за частотой выбранной гармоники.
Это если вам так хочется через Фурье делать. Слежение можно сделать и чисто временными способами. Зависит от ситуации. Если вы уже знакомы с Герцелем, то возможно через него вам будет проще.

Кроме того, т.к. sdft/герцель основаны на dft а не fft, вам никто не мешает выбрать нерегулярную частотную сетку, т.е. сделать собственный какой угодно шаг между частотными бинами, которые вы собрались анализировать. Возможно это упростит вам алгоритм интерполяции, хотя и не сделает его избыточным при заявленной точности.

Спасибо!

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

Я уже писал про параметры измерения частоты (25600 Гц, 320 мсек, окно Гаусса, FFT 8к), но вам такая точность, похоже, не нужна.
Вот результаты грубого моделирования с разной частотой самплинга на том же интервале 320 мсек Fft8k(25600Гц), Fft1k(3200Гц), Fft256(800Гц).
На входе синус 51 Гц + 10% третьей гармоники + 1% случайного шума.
1. FFT для трех случаев (8к, 1к и 256 точек)
Нажмите для просмотра прикрепленного файла
2. Обрезали все точки выше 75 Гц (чтобы не хватать вторую и третью гармоники, можно отбросить и точки ниже 25 Гц, но поленился это делать) и провели простейший фиттинг по гауссу. Показано для 8к.
Нажмите для просмотра прикрепленного файла

А вот численные результаты фиттинга. Обратите внимание на xc (положение центра, т.е. измеренная частота синуса) и его ошибку
Value Standard Error
A8k y0 0.23205 0.06367
A8k xc 50.99956 4.13929E-4
A8k w 9.05248 8.87164E-4
A8k A 28998.67448 2.76099
A8k sigma 4.52624 4.43582E-4
A8k FWHM 10.65848 0.00104
A8k Height 2555.94038 0.20737

Value Standard Error
A1k y0 1.03749 0.16077
A1k xc 50.99733 0.00105
A1k w 9.04926 0.00224
A1k A 28990.07709 6.9703
A1k sigma 4.52463 0.00112
A1k FWHM 10.65469 0.00264
A1k Height 2556.0913 0.52375

Value Standard Error
A256 y0 3.45667 0.4462
A256 xc 50.99655 0.0029
A256 w 9.04252 0.00622
A256 A 28920.49775 19.3375
A256 sigma 4.52126 0.00311
A256 FWHM 10.64675 0.00733
A256 Height 2551.85819 1.45431