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

[stepit]

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