Есть сигнал. Частота основной гармоники 50 гц.
Используется алгоритм ДПФ. При смещении частоты на 1-5 Гц. Появляется дополнительная погрешность по модулю ДПФ. Оконное сглаживание не годится.
частота дискретизации 1200гц задана жёстко и не меняется, выборок строго 24 за период.
Слышал про алгоритм интерполяции, но так ничего конкретного не нашёл.
Может кто сталкивался или известен источник.????
Необходимо минимизировать погрешность БПФ при смещении частоты сигнала при постоянной частоте дискретизации.

Посмотрите здесь
http://home.comcast.net/~kootsoop/EricJ2/index.htm
Там есть Матлаб-файлы.
Я когда-то пробовал их все. Macleod's estimator лучший в отношении систематических ошибок.
При достаточно большом отношении сигнал/шум они выходят все на предельные ошибки по Крамеру-Рао:
дисперсия ошибки измерения частоты CRLB(f) = 3/(2*pi*pi*N*N*N*SNR)
N - размер блока, SNR - сиигнал/шум
По мере увеличения размера блока и отношения сигнал/шум, когда предельная ошибка стремится к нулю, точность определяется систематическими ошибками алгоритма.

Да. такое случается.



Для чего? Чтоб померять мощность этой гармоники?



Погрешность? Типа значение на частоте 50 гц не соответствует амплитуде сигнала во времени?
Размер дпф какой?




Спорный тезис...




Чего интерполировать хотим? Спектр? Или его модуль?



Похоже, минимизировать тут нечего. Если частота гармоники не равна k*1200/N (N - размер дпф k = 0,1...N-1),
мощность гармоники будет размазана по 2-м ближайшим отсчетам дпф.

Будет размазано. Но если априорно известно, что в спектре присутствует только одна гармоника (на фоне шума) и известна инструментальная функция прибора - то достаточно определить с высокой точностью частоту (и соответственно энергию) в максимуме. Это возможно.
В десятки раз точнее, чем бин DFT. В соответствии с предельными оценками CRLB.
В простейшем случае достаточно вычислить "DFT" в точках отстоящих на пол-бина от энергетического максимума и провести интерполяцию параболой, с тем, чтобы найти максимум и его аргумент (частоту). В литературе это называется ML-extension. По приведённой ссылке другие более непрямые "оценщики"

Полностью согласен.

Интереснео, но не очень понятно. Похоже на передачу информации по силовой сети. Допустим, частота сети может быть 45-55 Гц. НО Она одинакова и для приемника и для передатчика. Так, что тут есть 24 точки для обоих. нормальная синхронизация и ничего не уходит и не размазывается.
"Вот уж, действительно, Все относительно, Все-все - ВСЕ!"

Это нужно для оценки самой сети (в трёхфазной сети скажем). нужно найти ДПФ модуль первой гармноники(50Гц). Нужно как-то обеспечить минимальное влияние плавания частоты на результат ДПФ.
Вообще изначально имеется сигнал 50Гц + гармоники (реальный сигнал сетевой частоты). В данном случае интересует информация об основной гармонике.

Один умный человек мне сказал:
Синхронизация хорошо осуществляться аналоговым способом. Строится
умножитель частоты на основе ФАПЧ. В вашем случае умножать надо на 24. Вот
эта частота с выхода умножителя и используется как частота квантования.
Эти импульсы запускают АЦП. В этом случае частота квантования становится
переменной.
С постоянной частотой квантования тоже можно, но если она выбрана с
огромным запасом. Например у нас есть примерно 360 точек в
периоде. Тогда мы можем подстраивать интервал по 1 кванту времени,
которому соответствует всего 1 градус. Такую подстройку можно и
программным путем сделать, подгоняя число используемых квантов времени под
интервал между переходами через ноль. Скажем, 358 или 361 вместо 360.
А с низкой и постоянной частой квантования... Не думаю... Путем
интерполяции, но это тоже связано будет со сложными вычислениями.
Вот про последнее хотелось узнать по-подробнее.

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

Так что, частота дискретизации задана от кварца, а частота сети "гуляет" и происходит набег фазы? Интересно. Если говорилось про интерполяцию, то это, возможно, для передискретизации. Только как ее делать, если соотношение частоты сети и 1200Hz может быть произвольным в любой момент. Теоретически, нужна PLL-ка на частоту сети, можно с удвадцатичетверенной частотой, только зачем тогда дискретизация строго на частоте 1200? Можно и на частоте PLL-ки. Если нельзя, то нужна передискретизация с применением интерполяции и децимации (+ фильтры) с переменными коэффициентами.

Если нужна амплитуда основной гармоники не нужно так извращаться, можно:
1) Считать энергетику сигнала не по одному бину FFT а по нескольким, т.е. оценить амплитуду на выходе полосового фильтра с полосой пропускания
[50Нz-delta ... 50Hz+delta]. В примитивном исполнении берете несколько квадратов модулей FFT вокруг интересующей частоты складываете извлекаете корень квадратный из суммы и получается оценка искомой амплитуды. Естественно это все работает если кроме интересующей гармоники в заданном диапазоне больше ничего не живет.

2) В топку FFT, делаете фапч , фапч захватывает искомую гармонику, считаете sqrt( avrg(sig*cos(w*t))^2 + avrg(sig*sin(w*t))^2 )
sig - входной сигнал, cos(w*t), sin(w*t) - выходы ГУН ФАПЧ и ГУН ФАПЧ сдвинутый на 90 градусов, avrg - усреднение сделанное тем или иным способом. Все это в цифре , естественно.
p.s. А вообще способ предложенный Fontp наверное самый эффективный...

Если частота плавает то уже не строго!

Этого не может быть интернет завален этим добром. Гугль на интерполяция выдает столько всего что полжизни не перечитать!

Что значит минимизировать погрешность БПФ? Если речь идет о размазывании гармоники, то оконное сглаживание помогает. Если точность оценки амплитуды не удовлетворяет при использовании оконного сглаживания, то тогда контур ФАПЧ можно предложить, тогда обработка с точностью до фазы даст требуемую точность оценки.

Что нужно оценить-то? Мощность первой гарионики? Дык накой тогда дпф нужно?
Нужно расчитать полосовой фильтр с f0=50 и полосой +-df где df - максимально допустимое отклонение частоты.
Можно даже квадратурный фильтр реализовать. И посчитать амплитуду/мощность сигнала на выходе этого фильтра.





И это чтоб померять мощность гармоники? Ну, тады ой...



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

... попробовал Macleod's estimator. Так и не понял точно что возвращается как x в результате вычисления и что делать дальше с полученным числом.
Это очень похоже на метод определения частоты по максимуму спектральной плотности...

Это он и есть. Только в десяток раз точнее, чем даёт максимум ДПФ, при разумном отношении cигнал/шум.
В широком диапазоне отношений сигнал/шум оценка выходит на предельно возможную по максимуму правдоподобия.
Если просто взять максимум спектра ДПФ, то максимальная погрешность определения частоты равна +- пол бина ДПФ.
Если найти максимум, а потом запустить тот Macleod's estimator, то он возвращает x-дробное смещение относительно максимального дискретного бина и как показывает модель точность в основном определяется шумом (и числом отсчетов N, конечно) и в разы лучше бина. То есть для единствиной гармоники можно победить дискретность сигнала и его ДПФ

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

Понятно, что всё это не работает, если линий в спектре 2 или несколько. Модель построена для одной спектральной линии на фоне белого шума.

как то это очень напоминает вот это ЭТО
Я лично не уверен, что статистические методы дадут результат на одном периоде сигнала.

Мы тут уже сделали прибор для измерения парметров сети на таком принципе. Только у нас частота оцифровки сильно выше (25600), но нам нужно было еще гармоники до 40-й считать. Так что, если только основную частоту - то может, и хватит. Точность алгоритма - порядка 1е-7, точность всего прибора (уже метрологическая) 0.01% и определяется, в основном температурным уходом резистров. Мы не используем скользящих буферов, берется буфер 8192 отсчета (отоло 16 периодов сети), на него накладывается Гауссовское окно (с правильной нормировкой, чтобы не изменяло выходных значений. Стандартная нормировка не годится), считается ДПФ, далее находится максимум пика, от него берется по несколько отсчетов в обе стороны (кол-во определяется параметрами Гаусса), суммирем их квадратично и получаем амплитуду гармоники. при таком алгоритме от частоты ничего не зависит в широких пределах. Единственная проблема - если частота изменяется на длине одного блока выборки - тогда ошибка увеличивается, но незначительно. В реальных измерениях этот эффект не приводит к большой погрешности по амплитуде. Вот с фазой - это да, там все начинает бегать очень сильно.

Если нужно что то померить и точности FFT не хватает, то можно порекомендовать такие варианты
1. Увеличить длину FFT, добив нулями входную последовательность - это увеличивает разрешение по частоте, выгодно использовать если нужно анализировать весь частотный диапазон
2. Использовать квадратичную и прочие интерполяции, уточняя положение максимума по трем бинам FFT - вроде так делают когда интересных точек в спектре не так много чтобы использовать [1], но и не так мало чтобы использовать [3]. В общем это некиий компромисс между экономией ресурсов и качеством, мне не очень нравится...
3. Поиск локального максимума спектра при помощи алгоритма Герцеля и любого из известных методов поиска экстремума. Способ хорош тем, что не накладывает априорных ограничений на точность, его удобно использовать при отсутствии шума (точнее когда шум сравним с машинной точностью), например если надо эксперементально проверить частотные характеристики рассчитанного фильтра. Если важна скорость, то для поиска лучше использовать метод Брента или другой с квадратичной интерполяцией. Понятно, что выгодно использовать если в спектре интересны 1-2 точки.

Все эти варианты увеличивают точность как по частоте, так и по амплитуде.

Если же интересует только амплитуда, то тут помогают окна без scalloping loss, типа FlatTop и прочих (Potter xxx, Mennen xxx). Ну а использовать их можно как с обычным FFT, если частота заранее неизвестна, или с тем же Герцелем. Последний вариант можно рассматривать как перестраиваемый фильтр - Герцель модулирует, а FlatTop - это тот же НЧ фильтр. Если частота известна с точностью до пол-бина то со стандартными окнами должно хватить одного прохода.

Добавление нулей НЕ УЛУЧШАЕТ разрешения спектрального анализа! Читайте матчасть прежде чем говорить. Иначе вы противоречите фундаментальному принципу неопредлнности Гейзенберга по малому промежутку можете получить сколь угодно высокое частотное разрешение!

Алгоритм герцеля ничем по сути не отличается от ДПФ, а следовательно вы не получите увеличения разрешения по частоте с его помощью. Алгоритм Герцеля применяют для рассчета спектральной составляющей на заданной частоте без необходимости вычисления ДПФ во всех точках. Принцип неопределенности он не не отменяет и разрешение не улучшает!


Очередное заблуждение! для увеличения точности оценки частоты необходимо увеличить интервал анализа во времени. Если не увеличивать интервал анализа то невозможно обеспечить улучшение разрешающей способности по частоте! Если у вас выборка сигнала во времени короче 1 секунды, то спектральные составляющие отстоящие менее чем на 1 Гц будут неразличимы что бы вы не делали.
Для точного измерения амплитуды используют сглаживающие окна.

Гейзенберг, герцы и секунды - это из физики, а вопрос судя по всему был про оцифрованный сигнал. Тут не то что Гейзенберга, а даже Котельникова уже поздно вспоминать, поскольку от сигнала осталась только последовательность x[i] из N цифр. И даже секунд между этими цифрами не осталось, а частота меряется уже не герцах или радианах в секунду, а просто в радианах (понятно о чем речь?).

А задача стоит как ловчее найти интересные точки на спектре, то есть на еденичной окружности z-образа последовательности x[i].
Заметьте, что хотя достаточно знать значения в N равноотстоящих точках на окружности (дискретный спектр), да и работать с ним конечно удобнее, никто на запрещает работать со значениями в произвольных точках (непрерывный спектр).

Например чтобы узнать значение z-образа в точке окружности с аргументом w можно в лоб посчитать сумму x[i]*exp(-I*w*i). Суммы такого вида называются полиномами, что легко увидеть, заменив exp(-I*w) на у. И считать полиномы можно разными способами, хоть в лоб, хоть Горнером, хоть Герцелем. Да, алгоритм Герцеля скорее ближе к правилу Горнера, чем к ДПФ Ведь ДПФ вычисляет значение полинома сразу на N равноотстоящих точках на единичной окружности, а Герцель считает значение в одной точке, зато произвольно выбранной. Только пожалуйста, не надо говорить что Герцель - это только то, чем DTMF декодируют. Ну да, декодируют, но ведь не только: погуглите Goertzel Horner ради интереса.

Ну так вот собственно о вопросе - ДПФ длины N дает значения не в любых точках спектра, а только в равноотстоящих N точках, включая единицу. На жаргоне эти точки называются бинами. Проблема в том, что интересные точки спектра (в данном случае речь идет о локальном максимуме вблизи заданной точки) могут и не совпадать с бинами... Хорошо еще если в соседних бинах значения спектра оказались равны - сразу можно догадаться что в силу симметрии максимум лежит посередине между ними.
А как быть если получатся например значения W[k-1] = 3; W[k] = 5; W[k+1] = 4?; Симметрии тут нет, ясно что максимиум между k и k+1, но вопрос, где именно? Можно просто добавить бинов, добив нулями вход (чтобы представить, как это работает, рассмотрите, как удлиненная последовательность умножается на матрицу ДПФ), но есть и другие подходы, которые собственно тут и обсуждаются.

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

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



ДПФ действительно определено так как Вы сказали. Но на самом деле частотно равноотстоящие гармоники можно было бы повернуть все вдруг на exp(iw0t) и ничего не изменится, эти N гармоник тоже будут базисом и ничего не случится )) А интересные точки попадут ровненько на получившиеся бины. Это тоже был бы спектр, но непривычный, не физический.
Поскольку мы не знаем заранее на сколько нужно повернуть )) то просто берут 3 точки вблизи максимума и проводят параболу. Получается хорошо.
Лучше бы было бы подогнать функцию окна под измеренные точки, но это слишком сложно вычислительно. (Если бы аккуратно подогнать SINC то получили бы ровно тот результат, который получается вставкой нулей - т.е. метод вставки нулей - это интерполяция синком. Образующиеся при вставке нулей новообразованые гармоники можно разложить по первоначальному базису и увидеть, что это именно так - интерполяция синком в частотной области).
Но поскольку верхушку функции окна в максимуме можно аппроксимировать параболой, то такой подход работает.

Вообще ДПФ лучше всего представлять себе как банк фильтров с равноотстоящими центральными частотами и одниковыми, но сдвинутыми частотными откликами. Тогда действительно отдельный фильтр из банка совершенно ничем не отличается от Герцеля.
Хотя принято ДПФ называть сразу все фильтры банка. Но нас в отношении измерения частоты изолированой гармоники интересуют только некоторые.

Лет 15 назад мы такое делали на 8051 микроконтроллере.
Меряли частоту через скользящий ДПФ с центром на 50 Гц.
И многие другие сигналы в рельсах метро.
Отклонение частоты от 50 Гц мерялось измерением периода сигнала с выхода скользящего ДПФ.
Если период =1 сек. - значит частота 49 или 51 Гц.
То или другое определялось по фазе.
Но никакой Герцель тут не катит из-за больших ошибок.
Амплитуда мерялась по максимуму интерполяции графика спектра в трех соседних точках.

А справедлив ли принцип неопределености для синуса?
А если сделать из конечной последовательности бесконечную путем периодического повторения, то увеличит ли это частотное разрешение?

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

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


Тем не менее, это все таки аппроксимация синк-интерполяции. Ее точность зависит как от окна, так и от метода, например quinn estimator работает с прямоугольным окном, но теряет точность для других.


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

Амплитуда мерялась по максимуму интерполяции графика спектра в трех соседних точках.

Объясните пожалуйста это делается...
Уважаемый fontp представлял коды матлаба(Macleods estimator), но там несколько непонятно...
И как можно искать частоту по максимуму спектральной плотности?

Типичная ошибка. По двум отсчётам (АЦП) невозможно (гарантированно) определить ни амплитуду, ни фазу, ни разность фаз самого наидеального синуса. Есть конечно предельный случай с вероятностью 0, когда эти два отсчёта совпадут с экстремумами синуса, но на это глупо уповать.

да по двум отсчётам это не вариант...

принцип неопределенности справедлив для всех процессов в известной нам вселенной.


НЕТ! Периодическое повторение синусоиды приводит к дискретизации спектра, но не улучшает разрешение по частоте.


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

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

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

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


Послушайте, уважаемый, вместо того чтобы устраивать бла бла шоу займитесь делом! Принцип неопределенности можно пояснить так: укорочение сигнала во времени НЕИЗБЕЖНО приведет к расширению спектра! Нельзя одновременно получить сколь угодно короткий сигнал со сколь угодно узким спектром! Если кто-то считает по-другому ну что ж... применительно же к квантовой физике, то принцип неопределенности играет такую же роль что и скорость света в теории относительности, и раскрытие неопределенности положение- импульс в квантовой механике можно сравнить с изобретением машины времени!

Принцип неопределённости, если не грузить его физическим или философским смыслом, всего лишь констатация свойства спектра Фурье.
С философской точки зрения можно сказать так: Принцип неопределённости имеет место всегда по отношению к разрешению, т.е. когда нужно различить две близкие спектральные линии. Но он не имеет никакого отношения к ситуации, когда априорно известно, что гармоника одна и нужно измерить её частоту. Измерять её мы будем в разы точнее, согласно статистическому критерию Крамера-Рао. В этом случае точность ограничена только отношением сигнал/шум.

Просто неверно точность измерения частоты одиночной гармоники называть разрешением. Давайте называть её погрешностью измерения, например



А как её ещё искать? ))

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

По ссылке, которую я привёл измеряют именно частоту комплексной экспоненты, причём с максимально возможной точностью. Там всё описано - измеряется спектр Фурье, ищется максимум спектра мощности, берутся 3 точки в окрестности спектрального максимума (n0) и формируется дробная поправка dn по 3-м спектральным отсчётам S(n0-1), S(n0), S(n0+1) с помощью интерполяционной формулы
Получают оценку частоты f = 2*pi* (n0+dn)/T

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

Здесь конкретно 50 Гц и фазовый шум - нестабильность турбогенератора - минимальный.
Сигнал до того сильный, что к нему помехи не прилипают, кроме как гармоники от нелинейностей.
Ну разве что там сварка, станки, лифты...
Если делать чистый скользящий ДПФ - то на выходе реальная и мнимая части - имеются натурально.
Скользящий ДПФ - та же узкополосная фильтрация, чем длиннее ДПФ - тем уже полоса.
А измерение 3-х гармоник в соседних ДПФ-каналах - это только для оценки амплитуды.
Опыт показывает, что такая оценка даже лучше, чем промышленным вольтметром переменного напряжения.

абсолютно верно! Чем длине ДПФ тем уже полоса и точнее оценка. Но изначально люди хотят на 20 мс получить оценку частоты с высокой точностью в интервале 45-55 Гц.

Это, очевидно, не то.
Действительная синусоида - это не одна спектральная линия, а пара. Причём, если они достаточно близко к нулевой частоте, то они начнут "сливаться" и рассматриваемые выше методы интерполяции работать не будут. А мы говорим о ПРЕДЕЛЬНО ТОЧНЫХ оценках, а не о каких попало.

Оценки, точность которых лучше чем по критерию Крамера-Рао невозможны, а которые заметно хуже нас не интересуют. Возникающие дополнительно проблемки при оценке частот реальных синусоид в отличие от комплексных (для которых вопрос давно раз и навсегда решён) рассматриваются отдельно до сих пор исследователями. Мне например попадалась такая статья, какая в аттачменте




Это не имеет отношения к делу. Мы обычно ищем максимально точную оценку частоты при возможно более коротком ДПФ

Дальше логика такая -

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

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

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



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



Мне кажется в статье "Optimum Windows for Carrier Frequency Estimation," http://lcs.syr.edu/faculty/sarkar/1983.html механизм влияния отрицательных частот изложен более понятно чем в прилепленной статье (см. параграф 'noiseless case' и fig.3). Наверное товарищи So и Chan эту статью не читали, иначе не стали бы занижать эффективность периодограммы, используя прямоугольное окно...



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

Принцип действительно могучий. Функция не может быть локализована и во временной и в частотной областях.
Т е если функция локализована во времени (дельта-функция), то спектр ее будет неограничен по частоте. И наоборот - локализация в частотной области даст неограниченный во времени сигнал.
Посчитать частоту гармонического дискретного вещественного сигнала x(n) без шумов c неизвестной амлитудой и фазой - не проблема. Достаточно 3-х последовательных отсчетов :

f=0.5*real( acos( (x(n)+x(n+2) )/(2*x(n+1))) )*Fd/pi


Проблема появляется при появлении шумов и помех.

Странный Вы человек, противоречивый.
So&Chan не пользуясь не нужными окнами и не читая ненужных статей вышли на теоретически предельно возможные по эффектиивности оценки частоты. Сама по себе периодограмма (без интерполяции) такую точность никогда не даст, какие окна не прикручивай. Ошибка будет оставаться порядка бина DFT. Окна всего лишь дают компромисс между шириной центрального лепества и боковиками.

А Вы говорите "занижать". Они не могли бы получить более точных оценок, чем по Крамеру-Рао, если бы даже на пупе извертелись ))
Хотя, конечно, могут быть разные подходы к интерполяции. И некоторые могут оказаться проще при реализации.

Странные какие-то споры.
Сигнал 50 Гц - это, в общем-то, стационарный сигнал, по крайней мере, на участке несколько минут.
Выбрать его период, пусть даже приблизителдьно, повторить сколько надо раз и анализируй вдоль и впоперек с какой хочешь точностью через ДПФ, БПФ с учетом его эргодичности.
Если период неизвестен, можно корреляцию посчитать с частотой дискретизации герц 500-2000.
Затем можно найти спектр по методу авторегрессии - как у Марпла написано.

Выкладываю свою тестовую среду для алгоритмов оценки частоты. Может пригодится тем кто кладет на принципы - пользуйтесь на здоровье, если разберетесь.
Сейчас там сравниваются метод поиска максимума периодограммы, несколько интерполяторов, и еще два алгоритма параметрических оценок, кстати из другой статьи того же So (параметрические естественно не будут работать если вход не соответствует модели 1 тон + шум).
Тест заключается в генерации тестового тона с заданными параметрами и с добавлением шума, затем частота тона оценивается алгоритмом и сравнивается с известной. Тест выполняется многократно для случайно изменяющейся частоты (в пределах 1-го бина) и начальной фазы, после усреднения получается оценка среднеквадратичной ошибки оценки частоты для заданного SNR. Такой тест выполняется для SNR от -6 до +inf dB, с использованием прямоугольного и других окон.



Нажмите для просмотра прикрепленного файла

Здесь таких нет. Поэтому кол-во скачиваний должно остаться равно нулю

Вы забыли написать про время измерения. 24 точки на период, а сколько всего точек?

... алгоритм выполняется постоянно, оценка БПФ производиться каждый период сигнала(20 мс), оценка может производится сколь угодно долго, главное - 24 точечное БПФ, которое должно выполниться за 20 мс, получается так... каждые 20 мс вычисляется БПФ... так вот, если частота начнёт гулять, нужно что бы данные полученные от БПФ (модуль) не были искажены...

Пытаться измерять частоту по одному периоду - это очень плохо. Тем более при помощи БПФ на 24 точки. Гармоники близко, они будут мешать интерполировать спектр.

Ну а хоть про то, что БПФ на 24 точки - это уже не Кули-Тьюки, а немного более сложная структура быстрого алгоритма, Вы в курсе?

А зачем в таком случае БПФ?
Тут наверное что-то следящее больше подходит, типа цифровой ФАПЧ.

Дык ведь, если доступен поток данных непрерывный. то и делаем интерполяционный фильтр с какой надо точностью (то есть, с необходимым нам увеличением частоты отсчетов, а потом анализируем. Проблем то нет. Лучше, все равно, не получится. Самое прямое и понятное решение. После интерполяции размер БПФ увеличенный берем, опять, такой какой нам нужен для точности.

Если надо получить оценку амплитуды каждые 20 мс делайте с перекрытием. Например БПФ на 64 точки, которое каждый раз сдвигается на 24 точки, т.е. постоянное перекрытие интервалов обработки. Таким образом длина БПФ отвязана от одного периода и можно делать весовое сглаживание для высокой точности амплитуды. Вот матлаб пример:

...цифровой ФАПЧ? нужно программное решение...



можно подробнее про этот подход...

На пальцах так.
Предположим (для упрощения) что достаточно увеличить число точек в 10 раз.
Делаем интерполяцию любым доступным инструментом (Matlab, Excel, ..).
Для этого добавляем по 9 нулей между каждым отсчетом и "получаем частоту дискретизации" 12000 Hz.
Пропускаем этот "поток" через цифровой фильтр нижних частот (интерполяционный фильтр), задача которого - обязательно пропускать первую гармонику (50 Hz), но достаточно хорошо () подавлять все частоты выше 600 Гц. Можно взять для пробы затухание 40-60 дБ, а там видно будет. Нужно брать FIR, так как он не трогает фазы.
Выходной сигнал будет иметь в 10 раз больше точек, что позволит повысить "точность".
Если есть шумы, а они есть, то применяемый фильтр их подавит (вне первой гармоники), поэтому чем длиннее фильтр, тем до некоторых пор это лучше. Нужно бы знать как быстро "плывет" сеть, чтобы выбрать фильтр, способный это отслеживать, если это надо.

Так это и есть программное решение.

... не совсем понимаю, что даёт интерполяция сама по себе. Пусть даже интерполировал не нулями, а вычисленными значениями (пусть интерполяция сплайнами).
Получил 24 выборки (анализ в течении 0.02с) -> Интерполировал, получил 240 точек на период -> провёл fft...
получил спектр, стало видно, что на погрешность это не влияет(((...
как стало видно увеличение числа точек даёт увеличение полосы спектра, но спектральные отчёты продолжают итди с тем же интервалом...

другое дело - это когда увеличиваем интервал анализа( скажем берём не один "период" (0,02с), а пусть 20 (0.4)), то тут спектральные составляющие идут уже с частотой 2.5 Гц (1/T), а не 50Гц.
Но проблема заключается в том, что нужно получать данные о БПФ каждые 20мс(стало быть по данным одного периода) и получается что априори имеем дело с 24 выборками((...

делайте обработку с перекрытием вычисляйте каждые 20 мс БПФ на интервале 0.4 секунды. Вот пример матлабовский: