Оконное взвешивание, алгоритм применения Опции

[Zelepuk]

Есть задача: нужно посчитать THD в промышленной сети 50Гц.
THD считается по 10 гармоникам (кратным оновной частоте). Основная частота меняется в пределах 45...55Гц.
Я хочу применить алгоритм:

-оцифровываю сигнал с частотой 4096Гц (это априорная величина задана жёстко),
-затем оконное взвешивание (использую flattopwin в MATLAB - окно с плоской вершиной),
-делаю БПФ на 512 точек (окно в 0,125с.)

Теперь самое интересное, если гармоник 10, то нужно просто найти 10 пиков в спектре и вокруг взять корень квадратный из суммы квадратов окружающих эти пики бинов. Но если частота основной гармоники плавает от 45 до 55 Гц, то получается например для 9-й гармоники бины на которых основная энергия меняются от 52 до 59 (при частоте основной гармоники 45-55Гц), а для 10-й Гармоники основная энергия распределяется на 58-70 бинах (в зависимости от изменения основной гармоники), 11-я гармоника уже бегает с 63 (при 45 гц основной гармоники) до 77 (при 55 гц основной гармоники) бина.

Иначе говоря трудно ожидать пики в одних и тех же местах спектра при изменении частоты основной гармоники, посему ищется способ однозначного определения местоположения пиков (амплитуд) гармоник в спектре. Речь идёт именно о поиске, так как на графиках точно видно что амплитуды всех гармоник имеют место быть с достаточно выскоой точностью.

[Alex11]

Позвольте не согласиться с т. bahurin по поводу того, что ничего не надо складывать. При наложении окна любая одиночная гармоника раскладывается на ряд бинов. Чтобы правильно подсчитать амплитуду исходной гармоники нужно взять корень квадратный из суммы квадратов бинов, на которые распалась гармоника. Далее нужно учесть коэффициент, даваемый окном, если только он не 1. Это зависит от нормировки оконной функции. Хотя в Вашем случае, когда нужно искать THD, т.е. отношение амплитуд, этот коэффициент сократится.
Полный алгоритм должен выглядеть следующим образом: ищем максимум основной гармоники в диапазоне 45-55 Гц, уточняем положение максимума по 3 или 5 точкам - ищем точное значение частоты (дробное). Вычисляем положения середин гармоник в бинах (целые, округленные). Вычисляем мощность основной гармоники как сумму квадратов бинов вокруг ее центра (сколько - зависит от окна, требуемой точности и уровня шумов). Вычисляем мощность остальных гармоник аналогично и складываем их. Затем делим одно на другое. Корень здесь брать не нужно, т.к. THD - отношение мощностей.
И еще одно замечание. Если используется комплексный FFT, то значение амплитуды бина - корень из суммы квадратов действительной и мнимой частей.

[Alexey Lukin]

Корень здесь брать не нужно, т.к. THD - отношение мощностей.

THD — это отношение амплитуд, а не мощностей. Корень брать нужно.

bahurin правильно написал: при использовании flat-top окна суммировать бины не надо, достаточно взять амплитуду максимального. При использовании других окон — надо суммировать.

Единственное что может быть плохо у flat-top окна в данном случае — это уровень боковых лепестков: он всего -70 дБ. Поскольку нелинейные искажения могут быть на уровне -100 дБ, flat-top окна может оказаться недостаточно!

[bahurin]

bahurin правильно написал: при использовании flat-top окна суммировать бины не надо, достаточно взять амплитуду максимального. При использовании других окон — надо суммировать.

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

Код

N = 512;

%время
t = 0:N-1;

%сигнал из 3 гармоник
s = 1.0*cos(2*pi*0.07*t)+0.5*cos(2*pi*0.153*t)+0.1*cos(2*pi*0.23*t);

%окно
w = blackmanharris(N);

%нормирую окно
w = N*w./sum(w);

%амплитудный спектр
S= abs(fft(s.*w))/N;

%частота
f = (0:N-1)/N;

% график
plot(f,S), grid;

Замените Блекмана с Харисом хоть флэттопом, хоть Хеммингом работать будет одинаково показывая амплитуду гармоники сигнала по максимальному локальному максимуму без всяких там суммирований и прочего шаманства.
Просто как бы я могу допустить, что можно суммированием вблизи гармоники сделать нормировку в частотной области аналогичную строке w = N*w./sum(w), но я не представляю как в частотной области сделать такую нормировку для сигнала со спектром не из одной гармоники а с более сложным спектром, например bpsk.

[Alexey Lukin]

Вот матлаб код который показывает спектр амплитуд сигнала без всякого суммирования как он есть по максимуму.
Замените Блекмана с Харисом хоть флэттопом, хоть Хеммингом работать будет одинаково показывая амплитуду гармоники сигнала по максимальному локальному максимуму без всяких там суммирований и прочего шаманства.

Ваш код работает плохо, неточно. Вся суть flat-top окна — в том, что с ним (и только с ним!) можно оценивать амплитуду гармоник по спектру без суммирования бинов с достаточной точностью (порядка 0.01 дБ). Для остальных окон вы получите ошибки от 0.7 до 4 дБ, т.к. амплитуда будет зависеть от частоты сигнала (см. хотя бы статью, которую выше рекомендует petrov).

Правильный метод оценивания амплитуды для других окон — суммирование энергии бинов в пределах каждого пика.

[bahurin]

Ваш код работает плохо, неточно. Вся суть flat-top окна — в том, что с ним (и только с ним!) можно оценивать амплитуду гармоник по спектру без суммирования бинов с достаточной точностью (порядка 0.01 дБ). Для остальных окон вы получите ошибки от 0.7 до 4 дБ, т.к. амплитуда будет зависеть от частоты сигнала (см. хотя бы статью, которую выше рекомендует petrov).

Правильный метод оценивания амплитуды для других окон — суммирование энергии бинов в пределах каждого пика.

1. Ну тогда просьба написать код который работает хорошо по вашему мнению. А так получается все суперэксперты а как доходит до подтверждения так лень 10 строчек в матлабе написать.
2. Как бы я не увидел ответа на вопрос что делать если сигнал не представлятся одной синусоидой, а имеет непрерывный по частоте спектр. В этом случае я не вижу возможности выделить отдельные бины (потому что их бесконечно много). Как тогда в этом случае оценивать спектр (если можно тоже с примером в матлабе плиз)?

Сообщение отредактировал bahurin - Aug 14 2011, 16:21

[Alexey Lukin]

1. Ну тогда просьба написать код который работает хорошо по вашему мнению. А так получается все суперэксперты а как доходит до подтверждения так лень 10 строчек в матлабе написать.

Я на Матлабе не шибко умею, может Петров вам подскажет? А что непонятно в моём объяснении?

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

Вообще непонятно, в чём вопрос. Если делается FFT, то спектр всегда дискретный. Если рассматривается бесконечный по времени сигнал, то спектр всегда непрерывный. От формы сигнала (синусоида/не синусоида) это не зависит.

Сообщение отредактировал Alexey Lukin - Aug 14 2011, 17:49

[bahurin]

Вообще непонятно, в чём вопрос. Если делается FFT, то спектр всегда дискретный. Если рассматривается бесконечный по времени сигнал, то спектр всегда непрерывный. От формы сигнала (синусоида/не синусоида) это не зависит.

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