Алгоритм расчета THD Опции

[stepit]

Здравствуйте, реализовал следующий алгоритм расчета 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, тогда пожалуйста поделитесь.

[SPACUM]

меня беспокоит точность определения амплитуды гармоники

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

[Alexey Lukin]

3.В Вашем случае невозможно точно синхронизировать частоту сигнала и частоту оцифровки, поэтому возможна только грубая оценка и если измерять THD+N не обязательно БПФ на много точек.

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

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

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

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

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

[stepit]

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

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

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

[Alexey Lukin]

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

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

- Почему, Вы, не используете chirp сигнал для расчета FR ? Ведь он должен дать лучшее разрешение по частоте и фазе.

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

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

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

Сообщение отредактировал Alexey Lukin - Aug 15 2011, 16:39