Здравствуйте, реализовал следующий алгоритм расчета THD:

1. Снимаем сигнал y со звуковой карты;
2. Накладываем окно Хенинга на сигнал: yw = hanning(y);
3. Делаем FFT и ищем частоту основной гармоники: S = fft(yw), f = max(abs(S));
4. Берем значения спектра в гармониках: h1 = S(f), h2=S(2*f).. hn = S(n*f);
5. Рассчитываем THD = 100 * norm([h2 .. hn]) / norm(h1).

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

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

Здесь про окна http://dsplib.ru/
здесь про алгоритм http://www.efo.ru/doc/Silabs/Silabs.pl?2089

1.Посмотрите как предоставляются соответствующие данные для реальных АЦП например для ADS1258. На графиках видно, что сигнал синхронен с оцифровкой и окон не используется, просто БПФ на 32к точек, где каждая гармоника представлена одной линией.
2.Для хороших АЦП шум и уровень гармонических искажений можно различить только после БПФ с таким большим числом точек. Интегральный уровень гармоник THD обычно меньше интегрального уровня шума. Регламентируется только THD+N.
3.В Вашем случае невозможно точно синхронизировать частоту сигнала и частоту оцифровки, поэтому возможна только грубая оценка и если измерять THD+N не обязательно БПФ на много точек.
4.Для грубой оценки получите логарифмический спектр, на нем виднее -> найдите выброс соответствующий поданной частоте и определите сумму квадратов амплитуд в него попавших -> разделите на нее сумму квадратов остальных амплитуд -> далее десятичный логарифм и умножить на 10. Это THD+N.
Можно и окна использовать, но все-равно грубо получится.
5.Для точной оценки нужна точная синхронизация, а если еще и много точек, то можно THD и N отдельно получить.

Спасибо за коментарии. Часто, в обзорах, встречаются графики зависимости THD от частоты. Вот и мне стало интересно как их расчитывают ? Делают несколько измерений на разных частотах, а потом интерполируют ? Или может посредством специальних сигналов (chirp, mls, multisin, etc.) ?
У кого какие мнения ?

1.Разделяю Ваш интерес, но кроме графиков в ПДФ на АЦП ничего интересного в литературе не попадалось.
Думаю есть установка делающая синус с THD+N < -150дБ и синхронно его семплирующая при помощи испытываемого устройства. Далее делается БПФ на достаточно большое число точек чтобы частоты THD были существенно больше среднего уровня шумов. (С увеличением числа точек шумовые составляющие уменьшаются пропорционально корню из числа точек, а составляющие THD не меняются.) Сделать можно, но не быстро.
2.Когда мне требовалось показать THD+N, я находил в выходном сигнале исходную синусоиду (по минимуму среднеквадратического отклонения) и вычитал ее. Вид оставшегося сигнала и спектр его как правило всех удовлетворяли.
3.Не верю в быстрые способы измерить THD на уровне -100дБ и сразу на всех частотах.

< Берем значения спектра в гармониках: h1 = S(f), h2=S(2*f).. hn = S(n*f);>
4.С окнами все будет точнее если учесть еще несколько частот Фурье около f, 2*f... (+- 3 или 4). Возвести амплитуды в квадрат просуммировать и извлечь корень.
А сколько получилось THD? А шума сколько?

Про зависимость THD от частоты в общем не слышал (он же по определению total!), но от частоты дискретизации - вполне вероятно, как и отношения с/ш. Есть такая программа SpecLab (от звуковухи работает), посмотрите на ней THD и с/ш при разных частотах дискретизации. При увеличении частоты дискретизации THD будет возрастать за счет добавления амплитуд шума на частотах "гармоник". Так что полосу анализа (и частоту дискретизации соответственно) надо выбирать ОЧЕНЬ осмысленно, чтобы не копить лишние шумы

Посмотрите, пожалуйста, работу Angelo Farina "Simultaneous measurement of impulse response and distortion with a swept-sine technique".



Тестировал алгоритм с помощью звуковой карты. Результаты сравнимые с Rightmark Audio Analyzer: THD = 0.0034%, THD+N = 0.0071%.

И еще, интересно, может кто знает opensource проекты где можно посмотреть как реализированы алгоритмы обработки сигналов: THD, THD+N, SNR, Clipping detection, etc.

Почему грубая? При достаточно большом размере БПФ гармоники в спектре чётко видны, и их амплитуды можно вычислить вполне точно (используя либо суммирование энергии вокруг каждой гармоники, либо flat-top окно). У хороших ЦАП и АЦП гармоники находятся на уровне -100...-110 дБ относительно уровня полезного сигнала, а спектр шума — на пару десятков дБ ниже (при размере БПФ порядка 500 мс).


Можно и интерполяцией, а можно и chirp-сигналом. Алгоритм тот же, что и для одной частоты. Учитываем только расширение спектра каждой гармоники из-за изменения частоты во времени.


Упомянутый вами RMAA — как раз и есть open source. Там открыты исходники измерительной части (правда, к старой версии).

Приблизительно так и делаю.
У меня прибор с АЦП ADS1254E, применял метод очень похожий, немного другое окно,
для измерения отдельно шума отрезал частоты гармоник +-6 гармоник Фурье и 6 гармоник Фурье в начале.
THD = -69дБ, N = -96дБ (для 80% максимальной синусоиды на 80Гц, генератор Г3-118).
Шум в начале спектра появляется при сигнале большем половины допустимой амплитуды.

1) Frequency response calculation. Как видно из кода, выходной сигнал у Вас состоит из набора синусоид разной амплитуды и частоты. А FR = db(InAmp) - db(OutAmp).
- За каким принципом, Вы, выбираете частоты и амплитуды для выходного сигнала (функ. FRGenerateOutput) ?
- Почему, Вы, не используете chirp сигнал для расчета FR ? Ведь он должен дать лучшее разрешение по частоте и фазе.

2) THD calculation.
- Какую weighting curve Вы используете при измерении THD+N (FFT::ACurve - не похоже что A-weighting curve) ?

Частоты — исходя из желаемой плотности точек графика и разделимости синусоид с помощью FFT фиксированного размера.
Амплитуды — согласно "среднестатистическому" спектру звуковых сигналов.


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


Для THD+N (A) используется именно A-взвешивание, согласно графику в вашей ссылке.