Я новерное не так вырозился.
У вас есть идеальная синусоида. Как вы на ней считаете гармоники их же там нет! Кроме основоной. Поэтому я и спросил что вы делаете с сигналом что бы появились гармоники.

в матлабе делаю так >S = sin(2*pi*t*53)+0.5*sin(2*pi*t*53*5)+0.3*sin(2*pi*t*53*3)+0.1*sin(2*pi*t*53*4);

это сигнал содержащий 5-ю, 3-ю и 4-ю гармоники.



n-ю гармонику относительно основной мы меряем когда надо THD померять. А вот если интересует именно амплитуда конкретной гармоники в сравнении с амплитудой этой гармоники в исходном сигнале?

У меня получается для сигнала S = sin(2*pi*t*53)+0.5*sin(2*pi*t*53*5)+0.3*sin(2*pi*t*53*3)+0.1*sin(2*pi*t*53*4); THD = 0.59 - это известно.
Когда вычисляю THD через fixed point FFT то получается 0,57.

Сообщение отредактировал Zelepuk - Aug 12 2011, 13:50

Сделайте printf от сигнала и проверьте, нет ли в нём клиппирования.

А потом проверьте допустимый входной диапазон для вашего FFT.

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

Зачем это делать в случае окна



?

Почему бы просто не найти 10 локальных максимумов в спектре и посчитать как отношение суммы квадратов амплитуд найденных 10-ти гармоник к квадрату амплитуды основного сигнала (для получения в процентах извлечь квадратный корень и умножить на 100)? Окно лучше использовать, конечно, иначе боковые лепестки спектра основного тона замаскируют гармоники. Возьмите окно с уровнем боковых лепестков -92 дБ (Блэкмана-Харриса, по-моему).

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

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

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

Есть флаттопы на любое подавление

http://www.rssd.esa.int/SP/LISAPATHFINDER/...ysis/GH_FFT.pdf

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



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

До тех пор, пока уровень гармоник на 33 дБ ниже уровня основного тона можно и Хэмминга использовать. Прилично иметь уровень боковых лепестков окна на ~10 дБ ниже минимального измеряемого сигнала. При оценке относительных значений компонентов спектра нормировка не требуется. Точность относительных измерений при большом отношении сигнал/шум определяется зависимостью амплитуды спектральной компоненты от частоты. При частоте сигнала, когда за анализируемый интервал времени поступает целое число периодов + половина, амплитуда в соседних бинах будет одинакова и минимальна. Если частота дискретизации в целое число раз выше частоты сигнала, оценка амплитуды будет максимальна. Окно с большим подавлением боковых лепестков уменьшает эту разницу в оценках амплитуды.

Сообщение отредактировал Дмитрий_Б - Aug 14 2011, 14:55

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

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

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

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

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


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

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

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

Спасибо, просветили!
Так в чём же вопрос заключается?