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

[TigerSHARC]

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

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

[fontp]

Есть сигнал. Частота основной гармоники 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 - сиигнал/шум
По мере увеличения размера блока и отношения сигнал/шум, когда предельная ошибка стремится к нулю, точность определяется систематическими ошибками алгоритма.

[thermit]

Есть сигнал. Частота основной гармоники 50 гц.

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

Используется алгоритм ДПФ.

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

При смещении частоты на 1-5 Гц. Появляется дополнительная погрешность по модулю ДПФ.

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

Оконное сглаживание не годится.

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

частота дискретизации 1200гц задана жёстко и не меняется, выборок строго 24 за период.
Слышал про алгоритм интерполяции, но так ничего конкретного не нашёл.

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

Может кто сталкивался или известен источник.????
Необходимо минимизировать погрешность БПФ при смещении частоты сигнала при постоянной частоте дискретизации.

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

[fontp]

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

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

[thermit]

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

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

[sup-sup]

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

Сообщение отредактировал sup-sup - Sep 4 2009, 18:25

[TigerSHARC]

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

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

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

Сообщение отредактировал TigerSHARC - Sep 5 2009, 12:30

[sup-sup]

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

[alex_os]

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

Если нужна амплитуда основной гармоники не нужно так извращаться, можно:
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 наверное самый эффективный...

[bahurin]

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

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

Слышал про алгоритм интерполяции, но так ничего конкретного не нашёл.

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

Необходимо минимизировать погрешность БПФ при смещении частоты сигнала при постоянной частоте дискретизации.

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

[thermit]

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

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

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

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

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

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

[TigerSHARC]

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

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

Сообщение отредактировал TigerSHARC - Oct 15 2009, 20:30

[fontp]

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

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

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

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

[bahurin]

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

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

[Alex11]

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

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

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

[анатолий]

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

Сообщение отредактировал анатолий - Oct 17 2009, 15:40

[AndeyP]

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

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

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

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

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

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

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

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

[TigerSHARC]

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

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

[GetSmart]

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

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

Сообщение отредактировал GetSmart - Oct 18 2009, 17:22

[TigerSHARC]

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

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

[bahurin]

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

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

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

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

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

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

[GetSmart]

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

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

[sup-sup]

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

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

[bahurin]

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

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

Принцип неопределённости придуман для охов

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

[fontp]

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

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

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

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

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

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

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

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

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

[анатолий]

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

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

[bahurin]

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

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

[fontp]

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

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

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

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

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

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

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

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

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

 

[sup-sup]

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

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

Сообщение отредактировал sup-sup - Oct 19 2009, 11:07

[AndeyP]

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

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

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

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

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

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

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

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

[thermit]

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

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

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


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

[fontp]

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

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

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

[анатолий]

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

[AndeyP]

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



 freq_est.zip ( 3.88 килобайт ) Кол-во скачиваний: 158

[GetSmart]

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

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

[Oldring]

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

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

[TigerSHARC]

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

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

[Oldring]

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

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

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

[alexkok]

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

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

[sup-sup]

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

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

Сообщение отредактировал sup-sup - Oct 22 2009, 04:14

[bahurin]

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

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

Код

clear all; clc;
fs = 1200;  %частота дискретизации
N = 1000;   % обработка в течении 1000 отсчетов
L = 64;     % окно анализа 64 отсчета
t = (0:N-1)/fs; %время
K = 24; %смещаем окно анализа на 24 отсчета
X =1;   % начальное смещение
i =1;   % номер итерации
f = 50+(rand(1,1)-0.5*10); %случайная частота 50+-5 Гц

% сигнал с небольшим гулянием амплитуды
s0 = 1.4*sin(2*pi*f*t).*(1+0.5*cos(0.5*pi*t+rand(1,1)*2*pi));
%цикл анализа
while((X+L)<N) % пока есть отсчеты для анализа
  
    s = s0(X:X+L-1);    %выбираю превое окно длины L = 64 отсчета
    w = window(@flattopwin ,L)'; % максимально плоское весовое окно чтобы получить максимально-точную оценку амплитуды
    S = abs(fft(s.*w))/(sum(w)); % БПФ  
    A(i) = max(S);  % амплитуда в данном окне такая вот
    T(i) = (X+L/2)/fs; %первое окно с задержкой 64 отсчета
    i =i+1; % инкремент номера цикла анализа
    X=X+K;% смещаю окно на 20 мс теперь будет без задержек
end;
plot(t,s0,T,A*2,'r.-'),grid; % синим сигнал, красным его амплитуда. точки оценки амплитуды идут через 20 мс.

[TigerSHARC]

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

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

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

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

Сообщение отредактировал TigerSHARC - Oct 27 2009, 15:07

[sup-sup]

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



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

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

[alexkok]

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

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

[TigerSHARC]

На пальцах так.
Предположим (для упрощения) что достаточно увеличить число точек в 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 выборками((...

Сообщение отредактировал TigerSHARC - Oct 28 2009, 17:41

[bahurin]

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

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

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

Код

clear all; clc;
fs = 1200;  %частота дискретизации
N = 1000;   % обработка в течении 1000 отсчетов
L = 64;     % окно анализа 64 отсчета
t = (0:N-1)/fs; %время
K = 24; %смещаем окно анализа на 24 отсчета
X =1;   % начальное смещение
i =1;   % номер итерации
f = 50+(rand(1,1)-0.5*10); %случайная частота 50+-5 Гц

% сигнал с небольшим гулянием амплитуды
s0 = 1.4*sin(2*pi*f*t).*(1+0.5*cos(0.5*pi*t+rand(1,1)*2*pi));
%цикл анализа
while((X+L)<N) % пока есть отсчеты для анализа
  
    s = s0(X:X+L-1);    %выбираю превое окно длины L = 64 отсчета
    w = window(@flattopwin ,L)'; % максимально плоское весовое окно чтобы получить максимально-точную оценку амплитуды
    S = abs(fft(s.*w))/(sum(w)); % БПФ  
    A(i) = max(S);  % амплитуда в данном окне такая вот
    T(i) = (X+L/2)/fs; %первое окно с задержкой 64 отсчета
    i =i+1; % инкремент номера цикла анализа
    X=X+K;% смещаю окно на 20 мс теперь будет без задержек
end;
plot(t,s0,T,A*2,'r.-'),grid; % синим сигнал, красным его амплитуда. точки оценки амплитуды идут через 20 мс.

[sup-sup]

... не совсем понимаю, что даёт интерполяция сама по себе ...
С помощью интерполяции (добавление нулей и фильтрация) мы получаем дополнительные точки для решения задачи во временной области. Очень важна фильтрация, которая позволяет убрать весь шум выше искомой частоты и точно расставить интерполяционные точки. Как вариант, фильтр должен делать свертку фрагмента выборки с добавленными нулями с импульсной характеристикой фильтра на каждый сдвиг, то есть с увеличенной Fs. Размер фильтра, как предлагал bahurin, может быть больше выборки. Допустим, коэффициент интерполяции равен 10, тогда точек в выборке будет 24х10=240. Размер фильтра может быть и больше, что желательно для лучшего подавления частот выше первой гармоники (то есть, можно подавлять со 100 Гц). Кстати, фильтр лучше сделать и до интерполяции, так как 24 точек хватит, тем более, что импульсную характеристику лучше взять длиннее. Еще раз, интерполяция нужна для решения задачи во временной области.
Насчет FFT на каждые 24 точки - это ведь не означает, что поток точек прерывистый. Применение FFT на 24 или на большее число точек нужно делать скользящим, т.е каждый раз со сдвигом на каждую точку (в 24 раза больше). После этого можно провести дополнительную обработку с целью получить пользу из этой избыточности.
TigerSHARC, Вы бы дали выборку для анализа, было бы легче. Только размер нужен, естественно, хотя бы в несколько периодов подряд.

Сообщение отредактировал sup-sup - Oct 28 2009, 19:12

[TigerSHARC]

... не совсем понимаю, что даёт интерполяция сама по себе ...
С помощью интерполяции (добавление нулей и фильтрация) мы получаем дополнительные точки для решения задачи во временной области. Очень важна фильтрация, которая позволяет убрать весь шум выше искомой частоты и точно расставить интерполяционные точки. Как вариант, фильтр должен делать свертку фрагмента выборки с добавленными нулями с импульсной характеристикой фильтра на каждый сдвиг, то есть с увеличенной Fs. Размер фильтра, как предлагал bahurin, может быть больше выборки. Допустим, коэффициент интерполяции равен 10, тогда точек в выборке будет 24х10=240. Размер фильтра может быть и больше, что желательно для лучшего подавления частот выше первой гармоники (то есть, можно подавлять со 100 Гц). Кстати, фильтр лучше сделать и до интерполяции, так как 24 точек хватит, тем более, что импульсную характеристику лучше взять длиннее. Еще раз, интерполяция нужна для решения задачи во временной области.
Насчет FFT на каждые 24 точки - это ведь не означает, что поток точек прерывистый. Применение FFT на 24 или на большее число точек нужно делать скользящим, т.е каждый раз со сдвигом на каждую точку (в 24 раза больше). После этого можно провести дополнительную обработку с целью получить пользу из этой избыточности.
TigerSHARC, Вы бы дали выборку для анализа, было бы легче. Только размер нужен, естественно, хотя бы в несколько периодов подряд.

... зачем нулями-то? Лучше же приближённо реальные точки найти... линейно, или сплайнами...
смоделирую в Матлаба завтра всё то, что вы сказали...
а насчёт выборки... ну вот примерно так в матлаб и задавать нужно

f = 50+(rand(1,1)-0.5*10); %случайная частота 50+-5 Гц
s0 = 1.4*sin(2*pi*f*t).*(1+0.5*cos(0.5*pi*t+rand(1,1)*2*pi));% сигнал с небольшим гулянием амплитуды

Сообщение отредактировал TigerSHARC - Oct 28 2009, 20:37

[sup-sup]

... зачем нулями-то? Лучше же приближённо реальные точки найти... линейно, или сплайнами...
смоделирую в Матлаба завтра всё то, что вы сказали...
а насчёт выборки... ну вот примерно так в матлаб и задавать нужно

f = 50+(rand(1,1)-0.5*10); %случайная частота 50+-5 Гц
s0 = 1.4*sin(2*pi*f*t).*(1+0.5*cos(0.5*pi*t+rand(1,1)*2*pi));% сигнал с небольшим гулянием амплитуды

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

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

[bahurin]

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

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

Вы только что сами все сказали интерполяция просто клонирует спектр. И если по исходному спектру нельзя точно померить, то и по "клонированному" тоже не получится, так как все повторяющиеся куски полностью совпадают с исходным. После фильтрации в идеале опять останется исходный спектр, по которому нельзя точно померить. А поскольку идеального фильтра нет, то вы только уменьшите точность оценки амплитуды. Еще раз повторяю: интерполяция не вносит никакой дополнительной информации, поэтому не может улучшить точности оценки.

Что касается фильтрации до интерполяции, то опять таки отсутствие идеального фильтра не полностью подавит высшие гармоники и использовать геометрические "точные" методы также не получиться. Этот вопрос уже обсуждался на форуме экспоненты.

[sup-sup]

Вы только что сами все сказали интерполяция просто клонирует спектр. И если по исходному спектру нельзя точно померить, то и по "клонированному" тоже не получится, так как все повторяющиеся куски полностью совпадают с исходным. После фильтрации в идеале опять останется исходный спектр, по которому нельзя точно померить. А поскольку идеального фильтра нет, то вы только уменьшите точность оценки амплитуды. Еще раз повторяю: интерполяция не вносит никакой дополнительной информации, поэтому не может улучшить точности оценки.

Что касается фильтрации до интерполяции, то опять таки отсутствие идеального фильтра не полностью подавит высшие гармоники и использовать геометрические "точные" методы также не получиться. Этот вопрос уже обсуждался на форуме экспоненты.

Так я ж не спорю, что интерполяция не генерирует дополнительную информацию - ей неоткуда взяться на том же числе точек. Интерполяция позволяет с максимальной точностью (используя предоставленную выборку) решить задачу на временной области (не на частотной, а на временной. И только). А предварительная фильтрация все равно улучшит сигнал (насколько сможет). Допустим, гармоника имеют -40 дБ (что хорошо), а фильтр -60. Наверняка хватит и отличаться от максимально возможного почти не будет.
На частотной области решения должны быть с той же точностью при условии использования того же временного окна. Для увеличения разрешения по частоте нужно выборку дополнить нулями и взять больший размер FFT (это тоже обсуждалось).
Если сравнивать FFT или свертку, то разницы конечного результата не будет при всех равных условиях.

[TigerSHARC]

Так я ж не спорю, что интерполяция не генерирует дополнительную информацию - ей неоткуда взяться на том же числе точек. Интерполяция позволяет с максимальной точностью (используя предоставленную выборку) решить задачу на временной области (не на частотной, а на временной. И только). А предварительная фильтрация все равно улучшит сигнал (насколько сможет). Допустим, гармоника имеют -40 дБ (что хорошо), а фильтр -60. Наверняка хватит и отличаться от максимально возможного почти не будет.
На частотной области решения должны быть с той же точностью при условии использования того же временного окна. Для увеличения разрешения по частоте нужно выборку дополнить нулями и взять больший размер FFT (это тоже обсуждалось).
Если сравнивать FFT или свертку, то разницы конечного результата не будет при всех равных условиях.

просто применил окно к сигналу, и получил погрешность чуть выше 1%

Код MATLAB:

T = 0.02;
Ts = 1/1200;
f = 55;
t = 0:Ts:T-Ts;
A = 10;
df = 1/T;
Fmax = 1/Ts;
n = length(Y);
w = hamming ( n, 'periodic');
f1 = 0:df:Fmax-df;
Y = A*sin(2*pi*f*t) ;
Xw = fft(Y.*w.')/sum(w);
X = fft(Y)/n;
subplot(2,1,1), stem(f1, abs(X)), grid
title('Спектр сигнала 24 выборки без окна')
subplot(2,1,2), stem (f1, abs(Xw)), grid
title('Спектр сигнала 24 выборки c окном')

[TigerSHARC]

... как ещё можно минимизировать погрешность???

[fontp]

... как ещё можно минимизировать погрешность???

Ещё можно адаптивной фильтрацией.
Вот здесь минимизируют через LMS предсказание синуса вместе с тремя гармониками
http://rapidshare.com/files/299646763/Nove...acking.pdf.html
Но следящая система будет работать не спешно

просто применил окно к сигналу, и получил погрешность чуть выше 1%

Y = A*sin(2*pi*f*t) ;
Xw = fft(Y.*w.')/sum(w);
X = fft(Y)/n;

Так это без шума. Без шума интерполяция в спектральной области даст почти бесконечную точность.
Сделайте преобразование Гильберта. Добавить в сигнал в конце нулей (от n до 3*n штук ) по вкусу, сделайте fft, найдите максимум, по 3-м точкам в окрестности максимума проведите параболу, найдите максимум этой параболы, аргумент максимума - это и есть частота

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

[sup-sup]

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

Интерполяция во временной области дает то же самое. Только нужно после интерполяции продолжать работать в той области, где она сделана. Кому что больше нравится (мне, как железячнику, более привычна временная).

[AlexU]

... как ещё можно минимизировать погрешность???

На входе почти синусоида, сиг/шум~4, диапазон 10..90 Гц, погрешность оценки частоты +/-0.02 Гц по спектру за 200 точек (частота выборок 200 Гц, время измерения 1 с). Все утопталось в ATmega8, осталось под индикатор и обмен с PC. Алгоритм используется полгода в другом серийном устройстве.

[GetSmart]

На входе почти синусоида, сиг/шум~4, диапазон 10..90 Гц, погрешность оценки частоты +/-0.02 Гц по спектру за 200 точек (частота выборок 200 Гц, время измерения 1 с). Все утопталось в ATmega8, осталось под индикатор и обмен с PC. Алгоритм используется полгода в другом серийном устройстве.

Это же противоречит принципу неопределённости!!!
Щас придут умные дяди и по голове настучат

[EvgenyNik]

TigerSHARC, Вы определитесь - Вам вычисление частоты на каждом периоде нужно или, всё-таки, пересчёт действующего значения основной гармоники?
Подозреваю, что нужно именно действующее значение основной гармоники. Вот его и считайте по 24 точкам.
Что касается частоты. Поскольку речь идёт о промышленной частоте, то надо понимать, что энергосистема - штука ну очччень инертная и за один период частота не уйдёт не то что с 50 на 49 (или 51), но и даже в пределах требуемой по стандарту точности определения, а значит такого разрешения по времени и не требуется. Считайте частоту, как минимум, на 10 периодах для защит и быстрой автоматики и (см. выше у AlexU) - 0,5 - 1,0 секунду для индикации.
Все видели кадры с Саяно-Шушенской? Железная "юла" там немаленькая, правда? А уж какая тяжёлая... Так вот, в штатном режиме такая турбина каждые 20 мс свою скорость не менять не будет... а вот динамические искажения сигнала и броски фазы быть могут (локальная нагрузка коммутируется и меняет местный режим сети довольно шустро, но не частоту!) - их то и будете отлавливать весьма успешно, но для дела бестолково.
Кстати, какова разрядность АЦП и какая точность определения частоты требуется? 0.01 Гц?

[AlexU]

Это же противоречит принципу неопределённости!!!
Щас придут умные дяди и по голове настучат

На работе прибора это никак не отразится.
Пока разговор вертится вокруг ограничений FFT и цитатами из учебников будем спокойно его выпускать.

[fontp]

Это же противоречит принципу неопределённости!!!
Щас придут умные дяди и по голове настучат

"Умные дяди" как раз утверждают, что это предельная точность, соответствующая критерию Крамера-Рао
Стандартное отклонение = примерно Fd/(SNR*NS**(3/2))
Можете сами проверить, что он именно её получил.

При измерении частоты во всём представимом диапазоне частот (от 0 до Fd/2) точность оценки улучшить нельзя,
поскольку критерий Крамера-Рао - это оценка максимального правдоподобия.
Если диапазон частот узкий, наверное, точность можно улучшить предварительной фильтрацией и передискретизацией. Только предварительная фильтрация вносит ещё и задержку, что противоречит вообще-то требованию, что измерение за минимальное время.

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

[AlexU]

"Умные дяди" как раз утверждают, что это предельная точность, соответствующая критерию Крамера-Рао
Стандартное отклонение = примерно Fd/(SNR*NS**(3/2))
Можете сами проверить, что он именно её получил.

При измерении частоты во всём представимом диапазоне частот (от 0 до Fd/2) точность оценки улучшить нельзя,
поскольку критерий Крамера-Рао - это оценка максимального правдоподобия.
Если диапазон частот узкий, наверное, точность можно улучшить предварительной фильтрацией и передискретизацией. Только предварительная фильтрация вносит ещё и задержку, что противоречит вообще-то требованию, что измерение за минимальное время.

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

Все верно, +1
Кстати, на практике увеличение частоты выборки не приводит к заметному увеличению точности. (Хотя диапазон соответственно расширяется)

Сообщение отредактировал AlexU - Nov 2 2009, 10:15

[GetSmart]

На входе почти синусоида, сиг/шум~4, диапазон 10..90 Гц, погрешность оценки частоты +/-0.02 Гц по спектру за 200 точек (частота выборок 200 Гц, время измерения 1 с). Все утопталось в ATmega8, осталось под индикатор и обмен с PC. Алгоритм используется полгода в другом серийном устройстве.

Откуда сигнал и какой он природы?

[AlexU]

Откуда сигнал и какой он природы?

Датчик (чего-то там ) частотный. На выходе от 7 до 3000 Гц. От 7 до 100 Гц сигнал очень слабый много шума s/n до 1,5. Бывают кратковременные импульсные помехи большой амплитуды. Форма сигнала зависит от частоты (5мВ синусоида <40Гц)-(15мВ синусоида ограниченная сверху или снизу, плавает постоянная составляющая <120Гц) -(20..500мВ синусоида<500Гц)-(3В почти меандр). Измерять по спектру - самое то. Хотя я бы для начала сам датчик переделал.

[TigerSHARC]

TigerSHARC, Вы определитесь - Вам вычисление частоты на каждом периоде нужно или, всё-таки, пересчёт действующего значения основной гармоники?
Подозреваю, что нужно именно действующее значение основной гармоники. Вот его и считайте по 24 точкам.
Что касается частоты. Поскольку речь идёт о промышленной частоте, то надо понимать, что энергосистема - штука ну очччень инертная и за один период частота не уйдёт не то что с 50 на 49 (или 51), но и даже в пределах требуемой по стандарту точности определения, а значит такого разрешения по времени и не требуется. Считайте частоту, как минимум, на 10 периодах для защит и быстрой автоматики и (см. выше у AlexU) - 0,5 - 1,0 секунду для индикации.
Все видели кадры с Саяно-Шушенской? Железная "юла" там немаленькая, правда? А уж какая тяжёлая... Так вот, в штатном режиме такая турбина каждые 20 мс свою скорость не менять не будет... а вот динамические искажения сигнала и броски фазы быть могут (локальная нагрузка коммутируется и меняет местный режим сети довольно шустро, но не частоту!) - их то и будете отлавливать весьма успешно, но для дела бестолково.
Кстати, какова разрядность АЦП и какая точность определения частоты требуется? 0.01 Гц?

мне очень интересно узнать вот такую вещ:
Есть сигнал (50Гц + гармоники). Вычисляем с помощью БПФ значения гармоник кратных основной частоте(50Гц).
Частота дискретизации фиксирована, временное окно наблюдения тоже. Проблема в том, что как только основная частота не равно 50Гц, то значения гармоник искажаются (оно и понятно, т.к. в окне наблюдения умещается не полное число периодов).
Основная задача - это узнать максимально достоверную информацию о значении гармоник.
При этом приоритетными задачами являются: постоянная(аппаратная) частота дискретизации, и минимальное временное окно наблюдение (насколько это возможно).

[EvgenyNik]

Если аппаратно, то речь, как я понимаю идёт о ПЛИС...
Тогда я Вам советую реализовать цифровой узкополосый КИХ-фильтр фильтр для 50 Гц (Матлаб и fdatool в помощь), который порежет всё остальное. К очищенному сигналу применить простой, но достаточно эффективный метод определения периода по переходам через ноль. И от этого "плясать" дальше.
Естественно, учитывать надо будет не 1 период, а 10-20 и момент перехода через ноль интерполировать (для дискретизации 1200 - линейно). Будет полезно сверять однородность периодов (хотя бы примерно) внутри интервала, чтобы избежать погрешностей от динамических искажений фазы.

[TigerSHARC]

Если аппаратно, то речь, как я понимаю идёт о ПЛИС...
Тогда я Вам советую реализовать цифровой узкополосый КИХ-фильтр фильтр для 50 Гц (Матлаб и fdatool в помощь), который порежет всё остальное. К очищенному сигналу применить простой, но достаточно эффективный метод определения периода по переходам через ноль. И от этого "плясать" дальше.
Естественно, учитывать надо будет не 1 период, а 10-20 и момент перехода через ноль интерполировать (для дискретизации 1200 - линейно). Будет полезно сверять однородность периодов (хотя бы примерно) внутри интервала, чтобы избежать погрешностей от динамических искажений фазы.

предполагалось использовать DSP-процессор.
А если требуется узнать информацию о большем количестве гармоник(20-40 гамоник), что бы на их основе коэффицент несинусоидальности посчитать...

[EvgenyNik]

коэффицент несинусоидальности посчитать...

Устройство контроля качества электроэнергии?
Вообще, на это дело есть пара иностранных стандартов и наш ГОСТ, где ширина окна наблюдения и интервалы усреднения оговариваются. Ищите...

(20-40 гамоник)

Ну тогда 1200 Гц это не просто мало, а мало до нереализуемости. Приглядитесь к частотам дискретизации 6кГц-10кГц

[TigerSHARC]

Устройство контроля качества электроэнергии?
Вообще, на это дело есть пара иностранных стандартов и наш ГОСТ, где ширина окна наблюдения и интервалы усреднения оговариваются. Ищите...

Ну тогда 1200 Гц это не просто мало, а мало до нереализуемости. Приглядитесь к частотам дискретизации 6кГц-10кГц

)) понятно что мало... пусть даже для предотвращения наложения используется RC-цепь и передискретизация...
Как минимизировать погрешность для всех гармоник (пусть 20 гармоник) при том что основная частота сигнала 45-55Гц... (т.е. нужно как то предотвратить растекание спектра...

[EvgenyNik]

для предотвращения наложения используется RC-цепь и передискретизация...

8-10-ая гармоника это максимум при 1200 Гц.

основная частота сигнала 45-55Гц...

Да ладно Вам, см. ГОСТ 13109-97

5.6 Отклонение частоты
Отклонение частоты напряжения переменного тока в электрических сетях характеризуется показателем отклонения частоты, для которого установлены следующие нормы:
- нормально допустимое и предельно допустимое значения отклонения частоты равны ± 0,2 и ± 0,4 Гц соответственно.

Если частота вышла за эти пределы, то коэффициент несинусоидальности как-то малоинтересен, т.к. в системе начинают происходить совсем другие процессы.

[TigerSHARC]

8-10-ая гармоника это максимум при 1200 Гц.

Да ладно Вам, см. ГОСТ 13109-97

Если частота вышла за эти пределы, то коэффициент несинусоидальности как-то малоинтересен, т.к. в системе начинают происходить совсем другие процессы.

Да знаю я теорему котельникова!!!

ладно объясню по порядку:
При ДПФ есть два источника методической погрешности:
1) наложение частот ввиду несоблюдения теоремы отсчётов;
2) неполное число периодов во временном окне;

ТАК ВОТ МЕНЯ ИНТЕРЕСУЕТ 2 ПУНКТ!!!!
пусть условия теоремы котельникова соблюдены абсолютно(эта задача уже решена)

Я говорю, то что знаю точно!!!
Существуют промышленные приборы(видел лично).
В руководстве по эксплуатации сказано, что выдаются значения до 40-й гармоники с такой-то погрешностью при том, что частота сети 45-55Гц.

[EvgenyNik]

Все эти приборы работают с окном наблюдения намного превышающем период сети, а потом результаты от нескольких замеров ещё и усредняются. Значение частоты, кстати, усредняется за время 20 секунд(!).
Если на пальцах, то выбрал один из двух способов определения частоты:
а) по переходам через ноль (в тяжёлых случаях - сигнал чистить)
б) из самого же Фурье - по приращению фазы основной гармоники. Грубо говоря, преобразовали на тех же 20 мс, получили модуль1-фазу1, на следующих 20 мс опять модуль2-фазу2 и т.п. Если частота сети идеальна, то разница фаза2 - фаза1 будет равна нулю. Если же частота меньше 50Гц или больше, то появится эффект скольжения частоты и Вы будете наблюдать поворот фазы в "+" или в "-". При 45(55)Гц это скольжение будет видно на двух соседних интервалах 20мс, при отклонении частоты на 0.05Гц период наблюдения увеличивается.
Допустим, с частотой разобрались. Начинаем изменять период таймера-задатчика частоты дискретизации в соответствии с определённой нами частотой, натягивая сетку квантования времени на реальный сигнал. При этом продолжаем "измерять" частоту.
В итоге добиваемся необходимой кратности частоты дискретизации и сетевой и все прочие расчёты остаются нетронутыми - никаких пересчётов коэффициентов, переменной длины буфера, комплексного доворота и т.п.

[TigerSHARC]

Все эти приборы работают с окном наблюдения намного превышающем период сети, а потом результаты от нескольких замеров ещё и усредняются. Значение частоты, кстати, усредняется за время 20 секунд(!).
Если на пальцах, то выбрал один из двух способов определения частоты:
а) по переходам через ноль (в тяжёлых случаях - сигнал чистить)
б) из самого же Фурье - по приращению фазы основной гармоники. Грубо говоря, преобразовали на тех же 20 мс, получили модуль1-фазу1, на следующих 20 мс опять модуль2-фазу2 и т.п. Если частота сети идеальна, то разница фаза2 - фаза1 будет равна нулю. Если же частота меньше 50Гц или больше, то появится эффект скольжения частоты и Вы будете наблюдать поворот фазы в "+" или в "-". При 45(55)Гц это скольжение будет видно на двух соседних интервалах 20мс, при отклонении частоты на 0.05Гц период наблюдения увеличивается.
Допустим, с частотой разобрались. Начинаем изменять период таймера-задатчика частоты дискретизации в соответствии с определённой нами частотой, натягивая сетку квантования времени на реальный сигнал. При этом продолжаем "измерять" частоту.
В итоге добиваемся необходимой кратности частоты дискретизации и сетевой и все прочие расчёты остаются нетронутыми - никаких пересчётов коэффициентов, переменной длины буфера, комплексного доворота и т.п.

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

[EvgenyNik]

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

[TigerSHARC]

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

Спасибо за совет! Думаю это то что надо!

[анатолий]

Если на пальцах, то выбрал один из двух способов определения частоты:
а) по переходам через ноль (в тяжёлых случаях - сигнал чистить)
б) из самого же Фурье - по приращению фазы основной гармоники. Грубо говоря, преобразовали на тех же 20 мс, получили модуль1-фазу1, на следующих 20 мс опять модуль2-фазу2 и т.п. Если частота сети идеальна, то разница фаза2 - фаза1 будет равна нулю. Если же частота меньше 50Гц или больше, то появится эффект скольжения частоты и Вы будете наблюдать поворот фазы в "+" или в "-". При 45(55)Гц это скольжение будет видно на двух соседних интервалах 20мс, при отклонении частоты на 0.05Гц период наблюдения увеличивается.
Допустим, с частотой разобрались. Начинаем изменять период таймера-задатчика частоты дискретизации в соответствии с определённой нами частотой, натягивая сетку квантования времени на реальный сигнал. При этом продолжаем "измерять" частоту.

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

[TigerSHARC]

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

Правильное решение - именно интерполировать сигнал относительно расчитанного значения основной гармоники.

[Den]

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

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

С уважением, Den.

[fontp]

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

С уважением, Den.

"достать" можно здесь
http://electronix.ru/forum/index.php?showt...mp;#entry726504

[Den]

Спасибо.

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

Сетевое напряжение 50 Гц. Необходимо вычислять фазу, частоту, амплитуду. Также необходимо проводить анализ гармонического состава.
Значение фазы, частоты и амплитуды необходимо получать как можно быстрее…

Планирую так:
АЦП 16 бит, частота дискретизации 6400Гц, FIFO буфер БПФ на 512 точек, т.е. получается 128 точек на период.
При поступлении каждой новой точки с АЦП в FIFO буфер, на него накладывается окно (Гауссовское?), производится БПФ. Далее нахожу амплитуду и фазу основной гармоники.
Остается частота основной гармоники, необходимая точность измерения которой равна 0,01Гц. Возможно ли её определение по скорости изменения мгновенной фазы основной гармоники? Или же необходимо построение параболы по максимумам?

С уважением, Den.

[TigerSHARC]

Спасибо.

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

Сетевое напряжение 50 Гц. Необходимо вычислять фазу, частоту, амплитуду. Также необходимо проводить анализ гармонического состава.
Значение фазы, частоты и амплитуды необходимо получать как можно быстрее…

Планирую так:
АЦП 16 бит, частота дискретизации 6400Гц, FIFO буфер БПФ на 512 точек, т.е. получается 128 точек на период.
При поступлении каждой новой точки с АЦП в FIFO буфер, на него накладывается окно (Гауссовское?), производится БПФ. Далее нахожу амплитуду и фазу основной гармоники.
Остается частота основной гармоники, необходимая точность измерения которой равна 0,01Гц. Возможно ли её определение по скорости изменения мгновенной фазы основной гармоники? Или же необходимо построение параболы по максимумам?

С уважением, Den.

получается вы хотите считать 128 точечное БПФ за 156,25 мкс. успеете?

Сообщение отредактировал TigerSHARC - Nov 21 2012, 07:43

[Den]

128 точечное БПФ у меня высчитывается за 1.1 мкс.
БПФ хочу делать по 4-м периодам, соответственно 512 точечное.


С уважением, Den.

[fontp]

Спасибо.

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

Сетевое напряжение 50 Гц. Необходимо вычислять фазу, частоту, амплитуду. Также необходимо проводить анализ гармонического состава.
Значение фазы, частоты и амплитуды необходимо получать как можно быстрее…

Планирую так:
АЦП 16 бит, частота дискретизации 6400Гц, FIFO буфер БПФ на 512 точек, т.е. получается 128 точек на период.
При поступлении каждой новой точки с АЦП в FIFO буфер, на него накладывается окно (Гауссовское?), производится БПФ. Далее нахожу амплитуду и фазу основной гармоники.
Остается частота основной гармоники, необходимая точность измерения которой равна 0,01Гц. Возможно ли её определение по скорости изменения мгновенной фазы основной гармоники? Или же необходимо построение параболы по максимумам?

С уважением, Den.

Амплитуда и фаза максимального бина будут скакать из-за размывания спектра.
Наоборот, сначала вы строите параболу по максимальному и соседним бинам. Определяете положение максимума (и возможно амплитуду).Потом уже гармонику точной частоты накладываете на входной сигнал и точно определяете фазу и амплитуду.
Известно, что такой способ определения частоты и фазы достигает предельно возможной точности (максимального правдоподобия), при заданной длительности сигнала и умереннных отношениях сигнал шум (0-30 дб). Большинство других методов, тот же ФАПЧ, не оптимальны с точки зрения измерений по коротким последовательностям. То есть точность может быть достигнута, но время захвата не соответствует критериям Крамера-Рао

Я приводил ссылки на серию исчерпывающих исследований Стэнфордского университета по поводу частного способа интерполяции спектра
с добавлением нулей
http://electronix.ru/forum/index.php?showt...amp;hl=Stanford

Там есть всё - выбор окон, влияние нестабильности частоты сигнала, и т.д. и т.п.

[Den]

Амплитуда и фаза максимального бина будут скакать из-за размывания спектра.
Наоборот, сначала вы строите параболу по максимальному и соседним бинам. Определяете положение максимума (и возможно амплитуду).Потом уже гармонику точной частоты накладываете на входной сигнал и точно определяете фазу и амплитуду.
Известно, что такой способ определения частоты и фазы достигает предельно возможной точности (максимального правдоподобия), при заданной длительности сигнала и умереннных отношениях сигнал шум (0-30 дб). Большинство других методов, тот же ФАПЧ, не оптимальны с точки зрения измерений по коротким последовательностям. То есть точность может быть достигнута, но время захвата не соответствует критериям Крамера-Рао

А разве наложение весового окна не устранит (не минимизирует) растекание спектра?

[fontp]

А разве наложение весового окна не устранит (не минимизирует) растекание спектра?

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

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

[Den]

Я приводил ссылки на серию исчерпывающих исследований Стэнфордского университета по поводу частного способа интерполяции спектра
с добавлением нулей
http://electronix.ru/forum/index.php?showt...amp;hl=Stanford

Там есть всё - выбор окон, влияние нестабильности частоты сигнала, и т.д. и т.п.

Большое спасибо, буду разбираться.

С уважением, Den.

Определяете положение максимума (и возможно амплитуду).Потом уже гармонику точной частоты накладываете на входной сигнал и точно определяете фазу и амплитуду.

Ещё один вопрос.

Под наложением гармоники точной частоты, Вы имеете ввиду свертку с комплексной синусоидой, частоты определенной по максимуму параболы?

[fontp]

Большое спасибо, буду разбираться.
С уважением, Den.
Ещё один вопрос.
Под наложением гармоники точной частоты, Вы имеете ввиду свертку с комплексной синусоидой, частоты определенной по максимуму параболы?

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

[Den]

Спасибо.
Буду пробывать.

С уважением, Den.

[Serg76]

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

[fontp]

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

Для одиночной комплексной синусоиды Rife и Burstin первыми дали статистическую оценку Крамера_Рао для одиночной синусоиды и
привели алгоритм, который в широких пределах параметров её достигает
Я загружал когда-то эту статью
http://electronix.ru/forum/index.php?showt...mp;#entry693650

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

Что значит "достаточно далеко"? Интуитивно можно догадаться, что это много дальше чем разрешение по критерию Релея 1/N
Насколько много не очень известно и зависит от многих деталей реализации.

Также существует теоретический результат, о том, что ситуацию нельзя в принципе улучшить, когда пики находятся достаточно близко.
Нижняя граница Крамера-Рао существует для любого многомерного случайного распределения. Это статистический критерий, который может , в принципе, быть получен в явном виде для любого многомерного нормального распределения. Такая граница была получена в том числе и для двух синусоид в белом гауссовском шуме. Результат очень громоздкий и не очень полезный для практических применений. Но ассимптотики он даёт такие как надо интуиции. Когда спектральные пики находятся далеко результат не отличается от критерия Крамера-Рао для одной синусоиды.
Когда пики приближаются к друг другу на расстояние соизмеримое с "разрешением" 1/N граница резко падает от
корень из(6/(N*(N-1)*(N-1)*(Es/No)))
к 1/N и хуже

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

непрактичный пруф для 2-х синусоид

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

 2_sinusoids_CRLB_best_worst.pdf ( 209.94 килобайт ) Кол-во скачиваний: 51

 

[Serg76]

2 fontp

Спасибо большое за такой развернутый ответ! интуитивно я так себе все это и представлял.
действительно, наверное с практической точки зрения проще не заморачиваться и использовать
однотональные алгоритмы, но при этом всегда надо помнить об "удаленности" паразитных гармоник
от основной. точность при этом, конечно падает, пробовал гонять разные однотональные алгоритмы
на случай полигармонического воздействия, во всех алгоритмах точность упала где-то на порядок в
рабочем диапазоне С/Ш при использовании прямоугольного окна. при других окнах точность еще
хуже, в отдельных алгоритмах уменьшение наблюдалось еще на два порядка. Но Маклеод лучший ,
как и рекомендуют буквари.

кстати, в рекомендуемой Вами работе Маклеода рассматривается работа "оценщиков" при
бигармоническом воздействии.

Еще раз большое спасибо.

[Serg76]

2 fontp

пока тема не "остыла", хотел еще спросить: помнится в одной из веток (не могу найти в какой ) Вы рекомендовали
при оценке частотной отстройки M-PSK сигналов использовать алгоритм Viterbi & Viterbi вместо возведения в степень.
Не могли бы Вы по-подробнее рассказать суть алгоритма и в чем его достоинство? Спасибо.

[fontp]

2 fontp

пока тема не "остыла", хотел еще спросить: помнится в одной из веток (не могу найти в какой ) Вы рекомендовали
при оценке частотной отстройки M-PSK сигналов использовать алгоритм Viterbi & Viterbi вместо возведения в степень.
Не могли бы Вы по-подробнее рассказать суть алгоритма и в чем его достоинство? Спасибо.

Смысл в том что фаза учетверяется, но амплитуда не возводится в степень. Вы гуглом можете найти много статей V&V
Я приводил эти.

http://electronix.ru/forum/index.php?showt...mp;#entry720790
Преимущество в том, что работает при более низком сигнал/шуме, но конечно слепая подстройка частоты конечно
очень сильно уступает подстройке по данным. Нужен очень длинный "пакет" данных, да и то ниже некоторого порога
сигнал/шум уже не накопляется даже увеличением длины пакета

2 fontp

Спасибо большое за такой развернутый ответ! интуитивно я так себе все это и представлял.
действительно, наверное с практической точки зрения проще не заморачиваться и использовать
однотональные алгоритмы, но при этом всегда надо помнить об "удаленности" паразитных гармоник
от основной. точность при этом, конечно падает, пробовал гонять разные однотональные алгоритмы
на случай полигармонического воздействия, во всех алгоритмах точность упала где-то на порядок в
рабочем диапазоне С/Ш при использовании прямоугольного окна. при других окнах точность еще
хуже, в отдельных алгоритмах уменьшение наблюдалось еще на два порядка. Но Маклеод лучший ,
как и рекомендуют буквари.

Лучший из простых оценщиков непосредственно по ДПФ. У любого алгоритма есть своя область применимости.
Маклеод устойчив к сильным шумам и у него из простых самая меньшая систематическая ошибка.
Но она есть, и она будет определяющей при высоком отношении сигнал шум > 40 дб

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

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

[Serg76]

Смысл в том что фаза учетверяется, но амплитуда не возводится в степень. Вы гуглом можете найти много статей V&V
Я приводил эти.

http://electronix.ru/forum/index.php?showt...mp;#entry720790
Преимущество в том, что работает при более низком сигнал/шуме, но конечно слепая подстройка частоты конечно
очень сильно уступает подстройке по данным. Нужен очень длинный "пакет" данных, да и то ниже некоторого порога
сигнал/шум уже не накопляется даже увеличением длины пакета

ок, спасибо. почитаю

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

Да, наблюдал это при моделировании.

Но она есть, и она будет определяющей при высоком отношении сигнал шум > 40 дб

к сожалению, о таком канале можно только мечтать

Кухня с интерполяцией спектра ДПФ в основном нашла применение для низкого сигнал-шума, в основном для модемов

мне это и надо Еще раз спасибо.

[Den]

Добрый день, уважаемый fontp!

Что-то я совсем запутался.
Накидал код в Matlab.


clear ;
%--------------------------------------------------------------------------
A = 50 ; % Amplitude 1-harmonics 50V , Идеальная амплитуда 1-й гармоники ;
Fmain = 50 ; % F 1 harmonics, Идеальная частота 1-й гармоники ;
%--------------------------------------------------------------------------
SNR = 80 ;
%--------------------------------------------------------------------------
Fs = 6400 ; % F sampling , частота дискретизации
N = 512; % количество отсчетов
%--------------------------------------------------------------------------
t = (0:N-1)/Fs; % вектор времени
freq = 0:Fs/N:200 ; % вектор частоты
%--------------------------------------------------------------------------
xx = A*sin(2*pi*Fmain*t) ; %
x = awgn(xx, SNR, 'measured');
%x = x + randn(size(x)) ;
%--------------------------------------------------------------------------
NFFT = abs(fft(x, N))./(N/2);
%--------------------------------------------------------------------------



subplot(2,1,1), plot(t,x), axis([0 t(N) -(1.4*A) (1.4*A)]) ;
hold on ;

subplot(2,1,2), stem(freq, NFFT(1:17)), grid on ;
hold on ;


subplot(2,1,1);
text(0, -20, sprintf('Идеальная частота %2.3f Гц\n',Fmain));

text(0, -50, sprintf('Идеальная амплитуда %2.3f В\n',A));



for k = 0.8:0.01:1.2
xx = A*sin(2*pi*(Fmain*k)*t) ;
x = awgn(xx, SNR, 'measured');
NFFT = abs(fft(x, N))./(N/2);

[Ampl_max N_bin_max] = max(NFFT) ; % максимальный бин и его значение

y_l = NFFT(N_bin_max - 1) ; % амплитуда левого бина
y_c = NFFT(N_bin_max) ; % амплитуда максимального центрального бина
y_r = NFFT(N_bin_max + 1) ; % амплитуда правого бина

x_l = (N_bin_max - 2)*(Fs/N) ; % частотное положение бинов
x_c = (N_bin_max - 1)*(Fs/N) ;
x_r = N_bin_max *(Fs/N) ;


inp_vector_x = [x_l x_c x_r] ;
inp_vector_y = [y_l y_c y_r] ;

out_vector_x = x_l:0.01:x_r ;

out_vector_y = lagrange(inp_vector_x, inp_vector_y, out_vector_x);


[ymax xmax] = max(out_vector_y) ;

A_calc = ymax ;

Freq_izmer = xmax ;



errorA = (abs(A - A_calc))*100/A ;


subplot(2,1,1) ;
h(1) = plot(t, x, 'r') ;
h(3) = text(0, -35, sprintf('Измереная частота %2.3f Гц\n', Freq_izmer));


h(4) = text(0, -65, sprintf('Измереная амплитуда %2.3f В\n', A_calc));
h(5) = text(0.009, -65, sprintf('( %2.3f )\n', errorA));






subplot(2,1,2);
h(2) = stem(freq, NFFT(1:17), 'r') ;
h(6) = plot(out_vector_x, out_vector_y, 'r') ;



drawnow ;
pause(0.5) ;
if k == 1.19
return
end
delete(h) ;
end




function yy=lagrange(x,y,xx)
% вычисление интерполяционного полинома в форме Лагранжа
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений аргумента, для которых надо вычислить значения полинома
% yy - массив значений полинома в точках xx

% узнаем число узлов интерполяции (N=n+1)
N=length(x);
% создаем нулевой массив значений интерполяционного полинома
yy=zeros(size(xx));
% в цикле считаем сумму по узлам
for k=1:N
% вычисляем произведения, т.е. функции Psi_k
t=ones(size(xx));
for j=[1:k-1, k+1:N]
t=t.*(xx-x(j))/(x(k)-x(j));
end
% накапливаем сумму
yy = yy + y(k)*t;
end



Получается, что построением параболы, амплитуду точно определить не получается?
Или я совсем уже туплю…
Ткните, пожалуйста, носом, в нужную сторону, а лучше конечно в пример Matlab.

С уважением, Den.

[fontp]

Добрый день, уважаемый fontp!

Что-то я совсем запутался.
Накидал код в Matlab.


clear ;
%--------------------------------------------------------------------------
A = 50 ; % Amplitude 1-harmonics 50V , Идеальная амплитуда 1-й гармоники ;
Fmain = 50 ; % F 1 harmonics, Идеальная частота 1-й гармоники ;
%--------------------------------------------------------------------------
SNR = 80 ;
%--------------------------------------------------------------------------
Fs = 6400 ; % F sampling , частота дискретизации
N = 512; % количество отсчетов
%--------------------------------------------------------------------------
t = (0:N-1)/Fs; % вектор времени
freq = 0:Fs/N:200 ; % вектор частоты
%--------------------------------------------------------------------------
xx = A*sin(2*pi*Fmain*t) ; %
x = awgn(xx, SNR, 'measured');
%x = x + randn(size(x)) ;
%--------------------------------------------------------------------------
NFFT = abs(fft(x, N))./(N/2);
%--------------------------------------------------------------------------



subplot(2,1,1), plot(t,x), axis([0 t(N) -(1.4*A) (1.4*A)]) ;
hold on ;

subplot(2,1,2), stem(freq, NFFT(1:17)), grid on ;
hold on ;


subplot(2,1,1);
text(0, -20, sprintf('Идеальная частота %2.3f Гц\n',Fmain));

text(0, -50, sprintf('Идеальная амплитуда %2.3f В\n',A));



for k = 0.8:0.01:1.2
xx = A*sin(2*pi*(Fmain*k)*t) ;
x = awgn(xx, SNR, 'measured');
NFFT = abs(fft(x, N))./(N/2);

[Ampl_max N_bin_max] = max(NFFT) ; % максимальный бин и его значение

y_l = NFFT(N_bin_max - 1) ; % амплитуда левого бина
y_c = NFFT(N_bin_max) ; % амплитуда максимального центрального бина
y_r = NFFT(N_bin_max + 1) ; % амплитуда правого бина

x_l = (N_bin_max - 2)*(Fs/N) ; % частотное положение бинов
x_c = (N_bin_max - 1)*(Fs/N) ;
x_r = N_bin_max *(Fs/N) ;


inp_vector_x = [x_l x_c x_r] ;
inp_vector_y = [y_l y_c y_r] ;

out_vector_x = x_l:0.01:x_r ;

out_vector_y = lagrange(inp_vector_x, inp_vector_y, out_vector_x);


[ymax xmax] = max(out_vector_y) ;

A_calc = ymax ;

Freq_izmer = xmax ;



errorA = (abs(A - A_calc))*100/A ;


subplot(2,1,1) ;
h(1) = plot(t, x, 'r') ;
h(3) = text(0, -35, sprintf('Измереная частота %2.3f Гц\n', Freq_izmer));


h(4) = text(0, -65, sprintf('Измереная амплитуда %2.3f В\n', A_calc));
h(5) = text(0.009, -65, sprintf('( %2.3f )\n', errorA));






subplot(2,1,2);
h(2) = stem(freq, NFFT(1:17), 'r') ;
h(6) = plot(out_vector_x, out_vector_y, 'r') ;



drawnow ;
pause(0.5) ;
if k == 1.19
return
end
delete(h) ;
end




function yy=lagrange(x,y,xx)
% вычисление интерполяционного полинома в форме Лагранжа
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений аргумента, для которых надо вычислить значения полинома
% yy - массив значений полинома в точках xx

% узнаем число узлов интерполяции (N=n+1)
N=length(x);
% создаем нулевой массив значений интерполяционного полинома
yy=zeros(size(xx));
% в цикле считаем сумму по узлам
for k=1:N
% вычисляем произведения, т.е. функции Psi_k
t=ones(size(xx));
for j=[1:k-1, k+1:N]
t=t.*(xx-x(j))/(x(k)-x(j));
end
% накапливаем сумму
yy = yy + y(k)*t;
end



Получается, что построением параболы, амплитуду точно определить не получается?
Или я совсем уже туплю…
Ткните, пожалуйста, носом, в нужную сторону, а лучше конечно в пример Matlab.

С уважением, Den.

Так Вы ничего и не считаете... Зачем Вам амплитуда центрального бина?
Но мы говорили о комплексной экспоненте. А не действительном синусе. От действительного синуса еще нужно как-то получить квадратуру. Такие статьи тоже приводились. Там
http://electronix.ru/forum/index.php?showt...66966&st=30
Для комплексной экспоненты берете любой метод из tst3.zip, лучше маклеод. Получаете интерполированое значение частоты.
Умножаете свой сигнал, на комплексную экспоненту полученой частоты комплексно сопряженную и суммируете для усреднения.
Или добавляете нулей , делаете дПФ, проводите параболу через три точки вокруг центрального бина. Берете не центральный бин, а верхушку параболы. Как здесь.
https://ccrma.stanford.edu/STANM/stanms/sta...14/stanm114.pdf

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

[Den]

Добрый день.

Пытаюсь всё разобраться. Вот ещё код.
Имею сигнал cod_ADC.


clear ;
%--------------------------------------------------------------------------
A = 50 ; % Amplitude 1-harmonics 50V , Идеальная амплитуда 1-й гармоники ;
Fmain = 51; % F 1 harmonics, Идеальная частота 1-й гармоники ;
%--------------------------------------------------------------------------
SNR = 90 ;
%--------------------------------------------------------------------------
Fs = 6400 ; % F sampling , частота дискретизации
N = 256 ; % количество отсчетов
%--------------------------------------------------------------------------
t = (0:N-1)/Fs; % вектор времени
freq = 0:Fs/N:200 ; % вектор частоты
[ee, N_points_freq] = size(freq) ;
%--------------------------------------------------------------------------

x_exp = A*exp(j*(2*pi*Fmain*t));
cod_ADC = A*cos(2*pi*Fmain*t) ;


%cod_ADC_complex = hilbert(cod_ADC) ;
cod_ADC_complex = cod_ADC ;


real_cod_ADC = real(cod_ADC_complex) ;
imag_cod_ADC = imag(cod_ADC_complex) ;


real_x_exp = real(x_exp) ;
imag_x_exp = imag(x_exp) ;


%--------------------------------------------------------------------------
subplot(3,1,1), plot(t,real_cod_ADC, 'r'), axis([0 t(N) -(1.4*A) (1.4*A)]), grid on ;
hold on ;
plot(t,imag_cod_ADC) ;
%--------------------------------------------------------------------------
subplot(3,1,2), plot(t,real_x_exp, 'r'), axis([0 t(N) -(1.4*A) (1.4*A)]), grid on ;
hold on ;
plot(t,imag_x_exp) ;
%--------------------------------------------------------------------------
NFFT_cod_ADC = fft(cod_ADC_complex) ; % Спектр
Amplituds_cod_ADC = abs(NFFT_cod_ADC)/(N) ; % Амплитуды

NFFT_x_exp = fft(x_exp) ; % Спектр
Amplituds_x_exp = abs(NFFT_x_exp)/(N) ; % Амплитуды

%--------------------------------------------------------------------------
subplot(3,1,3), stem(freq, Amplituds_cod_ADC(1:N_points_freq), 'r'), grid on ;
hold on ;

stem(freq, Amplituds_x_exp(1:N_points_freq)), grid on ;
hold on ;
%--------------------------------------------------------------------------
[Ampl_max_cod_ADC N_bin_max_ampl_cod_ADC] = max(Amplituds_cod_ADC) ;
[Ampl_max_x_exp N_bin_max_ampl_x_exp] = max(Amplituds_x_exp) ;

%--------------------------------------------------------------------------
piak3vect_cod_ADC(1:3)=NFFT_cod_ADC(N_bin_max_ampl_cod_ADC-1:N_bin_max_ampl_cod_ADC+1); % Isolated 3 samples around peak.
piak3vect_x_exp(1:3)=NFFT_x_exp(N_bin_max_ampl_x_exp-1:N_bin_max_ampl_x_exp+1);

%--------------------------------------------------------------------------
macleod_error_cod_ADC = macleod(piak3vect_cod_ADC) ;
macleod_error_x_exp = macleod(piak3vect_x_exp) ;

%--------------------------------------------------------------------------
% quin_error = quin(piak3vect) ;
% quadterp_error = quadterp(piak3vect) ;
% quin2_error = quin2(piak3vect) ;

min_er_x = Fmain - 0.1 ;
max_er_x = Fmain + 0.1 ;

figure ;


subplot(2,1,1), stem(Fmain, 1, 'r'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;
hold on ;
stem(((N_bin_max_ampl_cod_ADC-1)*(Fs/N))+(macleod_error_cod_ADC*(Fs/N)), 1, 'g'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;

subplot(2,1,2), stem(Fmain, 1, 'r'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;
hold on ;
stem(((N_bin_max_ampl_x_exp-1)*(Fs/N))+(macleod_error_x_exp*(Fs/N)), 1, 'g'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;


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

С уважением, Den.