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

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

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