минимизация погрешности, при плавающей частоте Опции

[TigerSHARC]

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

Сообщение отредактировал TigerSHARC - Sep 4 2009, 10:49

[AndeyP]

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

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

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

[bahurin]

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

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

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

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

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

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

[AndeyP]

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

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

Гейзенберг, герцы и секунды - это из физики, а вопрос судя по всему был про оцифрованный сигнал. Тут не то что Гейзенберга, а даже Котельникова уже поздно вспоминать, поскольку от сигнала осталась только последовательность 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, но вопрос, где именно? Можно просто добавить бинов, добив нулями вход (чтобы представить, как это работает, рассмотрите, как удлиненная последовательность умножается на матрицу ДПФ), но есть и другие подходы, которые собственно тут и обсуждаются.

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

[fontp]

Гейзенберг, герцы и секунды - это из физики, а вопрос судя по всему был про оцифрованный сигнал. Тут не то что Гейзенберга, а даже Котельникова уже поздно вспоминать, поскольку от сигнала осталась только последовательность 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 то получили бы ровно тот результат, который получается вставкой нулей - т.е. метод вставки нулей - это интерполяция синком. Образующиеся при вставке нулей новообразованые гармоники можно разложить по первоначальному базису и увидеть, что это именно так - интерполяция синком в частотной области).
Но поскольку верхушку функции окна в максимуме можно аппроксимировать параболой, то такой подход работает.

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

[AndeyP]

Здесь Вы не правы - принцип неопределённости связывает неопределённость локализации любой функции и неопределённость локализации её фурье-образа (забудьте Гайзенберга - вспомните оптику, дифракцию и определение разрешения).

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

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

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

Но поскольку верхушку функции окна в максимуме можно аппроксимировать параболой, то такой подход работает.

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

Лет 15 назад мы такое делали на 8051 микроконтроллере.
...
Но никакой Герцель тут не катит из-за больших ошибок.

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

[bahurin]

А справедлив ли принцип неопределености для синуса?

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

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

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

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

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

[GetSmart]

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

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

[fontp]

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

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

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

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

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

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

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

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

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

[анатолий]

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

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

[fontp]

Если делать чистый скользящий ДПФ - то на выходе реальная и мнимая части - имеются натурально.

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

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

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

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

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

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

Прикрепленные файлы

 Real_Freq_est.zip ( 689 килобайт ) Кол-во скачиваний: 317