Есть задача: нужно посчитать 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 гц основной гармоники) бина.

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

Гармоники в спектре ищутся как локальные максимумы: a[k-1] < a[k] >= a[k+1].
Выбираете несколько главных максимумов, превышающих по амплитуде уровень шума, — это и будут ваши гармоники.

А зачем вам flat top окно? Если вы считаете суммы квадратов, то оно только транжирит ваше частотное разрешение. Берите что-нибудь попроще и пошире, типа Blackman-Harris.

Я так и хотел делать - разбивать спектр на группы отсчётов и искать максимумы.
Дело как раз в том, что локальные максимумы находятся на разных расстояниях друг от друга. И при изменении частоты основной гармоники в широких пределах (45-55Гц) получается что группа для 8 гармоники (для 55Гц) "залезает" на группу для 9-й (для 45гц).

Остаётся только каким -то образом искать частоту основной гармоники и искать исходя из её значения. Или есть другой способ?


p.s. после применения окна Блэкмена-Харриса вижу, что точность определения гармоник ниже, чем прри использовании окна с плоской вершиной (flattop). Или нужно как-то просуммировать бины вблизи максимумов?

Вообще надо исходить из того с какой точностью вам надо померить уровни кратных гармоник и с каким динамическим диапазоном.

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

При этом я бы начал анализ частоты первой гармоники, потом предсказывал появление второй гармоники на удвоенной частоте (она точно не залезит на третью) потом предсказывал появление третьей гармоники и т.д. прошел бы столько сколько надо без проблем.

Я делал так для звука... Может поможет?

Не надо разбивать на группы. Просто ищите максимумы.


Вы в первом посте написали, что уже суммируете бины. Видимо, я неправильно понял.
Да, для других окон нужно просуммировать бины, принадлежащие каждому пику.

Можно сделать FIR фильтр в частотной области. Для этого нужен шаблон с основной частотой и гармониками (чем больше - тем точнее), который нужно 'сжимать - разжимать' в пределах диапазона частот 45-55 Гц. Размер и тип окна выбрать, чтобы только старшие частоты не сливались.

1. На группы нужно разбивать. Если частота первой основной гармоники может находится в пределах 45-55 Гц, то после взятия FFT необходимо искать локальный максимум в диапазоне 45-55 Гц. Никакие бины не складывать, просто методом перебора найти максимальный. Найдя этот локальный максимум надо запомнить его уровень A0 и его частоту f0. Для увеличения точности оценки максимума и частоты можно по трем точкам уточнить по полиному второй степени.

2. После того как мы оценили частоту основной гармоники, то предсказываем появление второй гармоники на частоте 2*f0 +- df, где df есть диапазон выбора второго локального максимума определяется точностью оценки частоты первой гармоники. Например если у вас частота дискретизации 4096 Гц, и вы берете 1024 точки, то точность оценки первой гармоники будет +-2 Гц. Тогда можно для второй гармоники взять df = 6 Гц (ошибка в 2 Гц также удвоится и будет +-4, с запасом можно взять +-6). Аналогично п.1 ищем локальный максимум, ничего не суммируя. Получили уровень второй гармоники равный A1 и частота второй гармоники f1. ВАЖНО частота второй гармоники также определена с точностью +- 2 Гц.

3. Предсказываем частоту третьей гармоники равную f2 = 3*f1/2. Поскольку мы вторую гармонику оценили с точностью +-2 Гц, то df для третьей гармоники можно также взять уже равной 4 Гц, поскольку df = 3*(+-2)/2 = +-3, а мы возьмем +-4, с запасом. Снова ищем локальный максимум запоминаем его уровень A3 и частоту f3

4. Для четвертой гармоники предсказываем частоту ее появления на f4 = 4*f3/3. далее везде полосу для поиска локального максимума можно брать +-4 Гц.

Итак ничего не надо суммировать и т.д., нужно лишь запусить итерационную процедуру поиска локальных максимумов в узкой полосе предсказывая каждый раз очередную гармонику. Соотвественно потом надо поделить полученные уровни An на A0.

Я эспериментирую в CBuilder над FFT и окнами.
Нашёл реализуцию fixed-point FFT http://www.jjj.de/fft/fix_fft.tar.gz

Генерирую в Матлаб:
- сигнал (чистая синусоида частотой 53 гц - чтобы увидеть растекание спектра) длиной 1024;
- окно с плоской вершиной длиной 1024;

Все массивы в матлаб перевожу в int16 и копирую в CBuilder.

Делаю комплексное фурье(с использование подпрограммы на ссылке выше). Нахожу амплитуду (корень квадратный из суммы квадратов мнимой и действительной части поэлементно соответственно).

В итоге получаю через printf нечто вроде

0
0
0
16
1256
3256
1245
145
11
0
1
0
1
0
0
0
1
0

(все числа вымышлены и только для наглядности).

Когда суммирую числа из кучки выше (16+1256+3256+1245+145+11) и умножаю на 2 - получается в точности то, что и в матлаб!
Даже удивился точности!
Значит всё работает подумал я!

Да не тут то было.

Как только к сигналу примешиваю 5-ю гармонику с амплитудой 0.5 от основной, то сразу же спектр, выдаваемый целочисленным алгоритмом в CBuilder сильно искажаться по всем гармоникам, в итоге получаю кучу помех в спектре, причём там где спектр должен быть нулевым.

например

0
10
1
1
16
32
1114
123
0
1
0
1
10
32
556
124
11
0
1
0
1
15
123
112
144
114
10
1
1



Но в матлаб всё работает превосходно.

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

Почему так может происходить?
Может дело в реализации FFT? Может кто поделится другой реализацией (только для int16). На плавучку очень не хочется переходить (код будет перенесйн на микроконтроллер без FPU)

меня сильно настораживает тот факт что вы так спокойно решили складывать и умножать подгоняя свои результаты под матлабовские. Миром правят нормировки. По счастливой случайности складывая "кучку" вы получили одну и туже нормировку с матлабом. В другом случае у вас не проканало. Может стоит сначала отмоделировать все в матлабе и потом пошагово искать и устранять разночтения в свое программе

к сожалению с целыми числами полноценно работать в матлабе не научился.

Можно подробнее про нормировку?

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

Zelepuk а как вы получаете гармоники? Из нулевой гармоники.

всмысле как? Я просто умножаю на окно, делаю БПФ и затем смотрю что получилось. Нулевая гармоника - это постоянная составляющая.



Дело как раз получается в том, что я делаю всегда одно и тоже с сигналом: сигнал+окно-> БПФ -> суммирование на некотором интервале (именно суммирую, я просто складываю "палки" и всё)

Если сигнал есть одна гармоника (пробовал разные частоты от 45 до 55Гц) - всё красиво и точность по амплитуде 0,01%, если сигнал представляет собой сумму гармоник (пусть есть 7-я и 5-я гармоники с амплитудами 0,1 от основной) - наблюдается погрешность.

По вашим словам получается когда сигнал чистая синусоида, то мой алгоритм учитывает ослабление от окна, когда сумма гармоник - не учитывает. Но я всегда делаю с сигналом одни и теже операции.

нет я говорю о том уже в котором посте, что не надо никакие гармоники складывать, если вы хотите померить уровень сигнала, надо мерить его по максимуму при учете к-та ослабления окна. Если вам надо знать уровень n-ой гармоники относительно основной, то и учитывать ослабление окна не требуется при взятии отношения оно само исчезнет.

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

в матлабе делаю так >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.

Сделайте 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 дБ ниже минимального измеряемого сигнала. При оценке относительных значений компонентов спектра нормировка не требуется. Точность относительных измерений при большом отношении сигнал/шум определяется зависимостью амплитуды спектральной компоненты от частоты. При частоте сигнала, когда за анализируемый интервал времени поступает целое число периодов + половина, амплитуда в соседних бинах будет одинакова и минимальна. Если частота дискретизации в целое число раз выше частоты сигнала, оценка амплитуды будет максимальна. Окно с большим подавлением боковых лепестков уменьшает эту разницу в оценках амплитуды.

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

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

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

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


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

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

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

Повторяю еще раз. Как производить оценку уровня спектральных составляющих, если нет возможности суммировать отдельные бины гармоник, потому что нет в сигнале отдельных гармоник, а есть непрерывный широкий спектр, продискретизированный по сетке частот на выходе fft? Поскольку алгоритм спектрального анализа должен работать независимо от того сколько гармоник в сигнале и как близко они друг к другу стоят, то совершенно не ясно как производить суммирование, когда они вообще друг с другом сливаются.

Если это сливающиеся синусоиды и нет возможности увеличить частотное разрешение анализа — надо применять параметрические методы.

Ну понеслась. Вы вообще себя слышите? Есть сигнал я не знаю что это за сигнал и не знаю вообще о нем ничего (я даже не знаю вообще есть ли он). Есть анализатор спектра. Я подаю свой сигнал на анализатор спектра (а может и не подаю, потому что забыл питание включить). И как он понимает что ему надо складывать гармоники, или не надо складывать, или он должен какие-то параметрические методы применить? И как по вашему в спектроанализаторе происходит оконное взвешивание? Я так понимаю дискуссия переходит в разговор слепого с глухим. Вы не ответили ни на один из моих вопросов, не привели ни одного внятного примера. Я умываю руки.

А я сделал в железе. Результат отвратительный, особенно для флаттоп. Если применяется любое окно кроме прямоугольного, то энергию входного сигнала и энергию спектра невозможно связать никакой формулой. То же и для максимумов пиков спектра. Даже для одной частоты.
Для синуса одна формула, для косинуса другая, а если еще частота не совпадает?...

Alex11 говорил о том, что сделано в железе с применением окна Гаусса, точность 0,01% по амплиуде. Нужна лишь длинная выборка.

Значит, где-то ошиблись. Должно быть соответствие как для общей энергии спектра, так и для максимума.

Шутите. Для синуса подавляются малоэнергетические края, а для косинуса самые максимумы. Для любых окон так. Не работаю на матлабе, все уже в СИ переписано.

Это что означает?

Пусть признается сколько гармоник и на какой основная частота. Проверить просто: синус, косинус и погрешность готова.


Возьмем один период синуса. Умножим на окно. Найдем среднеквадратичное значение.
Сделаем то же для косинуса. Будет обязательно меньше. И значительно.

Ну и что? На выходе максимального бина для флаттоп значение амплитуды будет одно и то же вне зависимости от частоты и фазы.

Проверяли?
У меня коэффициенты:
a0 = 0.21557895;
a1 = -0.41663158;
a2 = 0.277263158;
a3 = -0.083578947;
a4 = 0.006947368;
Может неверно?

Даже проверять нечего. Поймите простую вещь про FFT, можно абстрагироваться от быстрого алгоритма и представить его как набор N комплексных гетеродинов, за ними N ФНЧ, за ними N дециматоров в N раз. Без окна ФНЧ это просто скользящее среднее, для флаттоп собственно ФНЧ с ИХ флаттоп. Подаём синус, на выходе демодулятора с максимальной амплитудой в случае флаттоп разумеется ничего не зависит от частоты, поскольку вершинка фильтра плоская, и уж тем более от фазы ничего зависеть не может, ну повёрнут/вращается комплексный вектор, квадратуры меняются, а амплитуда постоянна.

Перепроверю.

SPACUM и Петров, вы оба по-своему правы.
Ошибка SPACUM в том, что он пытается измерить слишком низкую частоту для данного размера окна. Следует помнить, что flat-top (как и другие окна) даёт аккуратную оценку амплитуды только тогда, когда в пределы главного лепестка попадает один единственный тон. В примере SPACUM с синусом и косинусом в пределы главного лепестка очевидно попадает 2 тона: f и -f.
Вывод: для аккуратной оценки нужно, чтобы частота сигнала была не ниже, чем полуширина главного лепестка окна.

Кажется понял что имеется ввиду, при анализе действительного сигнала, когда в один бин попадают и положительная и отрицательная частоты, на выходе фильтра наблюдаются биения, тут конечно нужно разрешение увеличивать.

Был не прав. Прошу прощения. Амплитуды выбросов АЧХ окна флаттоп на синусах и косинусах равны.
При оцифровке 12кГц и числе точек 2000 на висящем куске провода на логарифмическом спектре раздельно видны 50 гармоник.
Нечетные стоят, четные меняются.
Идеи:
1. Продифференцировать сигнал до АЦП чтобы совсем мелких гармоник не было.
2. На логарифмическом спектре использовать кепстр для определения 1/ fundamental frequenсy

Коли спрашивали, признаюсь: оцифровка на частоте 25600 Гц, основная измеряемая гармоника - 50 Гц. При вычислении RMS используются все до Fs/2, при расчете гармонического состава - 40 гармоник. Но в последнем случае очень тяжело смотреть погрешности, т.к. все существующие эталоны дают отвратительную точность. То, что приводилось (0.01%) - это точность измерения RMS и амплитуды основной гармоники. Это на существующем приборе. При моделировании точность рассчетов (смесь float и double) на уровне 1e-8.

Большое спасибо за ответ. Мне такую задачу приходится решать при калибровке приборов и именно из-за очень плохих эталонов. В виброметрии других не бывает. Попробую Ваш метод. Но только в конце сентября. Работы много и отпуск. В реальном приборе точность расчетов (int32) на уровне 1e-8 для числа точек 2^n (на уровне 1.5e-7 для числа точек 2^n*10^n).

Здесь было предложено уже достаточно грамотных решений на основе преобразования Фурье. Я не буду повторять их. Также не буду говорить про затратные альтернативы типа MUSIC (мне самому стало интересно, как его упростить или построить новый алгоритм для случая кратных гармоник, может получится и не очень затратно). Позволю себе пофантазировать.
В описании сигнала ничего не сказано о шуме. Слукавлю и буду считать его несущественным (он таки может быть несущественным по сравнению с нелинейностью). Тогда достаточно сравнить мощность первой гармоники с мощностью остального сигнала без учета шума. Заметьте, знать точно частоту "основного тона" и мощности каждой гармоники в этом случае не нужно.

Потребуется три фильтра (не считая ФНЧ перед/после АЦП). Первый (простейший рекурсивный без умножителей) вырезает "постоянку" (смещение АЦП), не трогая первую гармонику и выше. Второй может быть вычислительно совмещен с третим (зависит от реализации), основное назначение - разделить сигнал на две полосы: ниже 55 и выше 90 Гц. Как вариант, при независимой реализации ФНЧ/ФВЧ, выход КИХ ФНЧ можно считать сразу децимированным (опять же для экономии). Далее усредняем мощность с выходов двух фильтров и делим. Либо, что экономнее, итеративно находим усреднение отношения двух мощностей без деления. Ценой одного умножения этим же методом можно получить и квадратный корень отношения двух мощностей. Либо, если хочется иметь значения в dB, то вычислить разность двух логарифмов. Возможно, потребуется согласовать задержки фильтров, чтобы скомпенсировать ошибки при изменении амплитуд. Естественно, точность измерения будет зависеть от избирательности фильтров (это верно и при выборе оконной функции при использовании преобразования Фурье), неравномерности в полосе пропускания, длины окна усреднения, точности арифметики. Ну и конечно же от шума, пренебрежение которым было разрешающим для подхода. При непрерывной (поотсчетной) обработке использование рекурсивных фильтров должно быть значительно экономнее по затратам, однако... однако не буду отвлекаться на вопросы реализации рекурсивных фильтров на малобитной арифметике (MCU в сопутствующей теме назван слабым).

Фантазирую дальше. Первый ФВЧ, удаляющий постоянную составляющую, можно совместить с последующим ФНЧ, что открывает возможность для добавления Гильбертовой пары для обоих полос, и от непрерывной обработки можно перейти к блочной. Для неперекрывающихся окон затраты, не считая накладных расходов, окажутся: 4 умножения и 4 сложения на отсчет для фильтров, 4 умножения, 2 сложения и одно деление (или два логарифма) на блок для вычисления мощностей и их отношения. При этом затраты памяти - одна ячейка для текущего отсчета и 4 аккумулятора по одному на каждый фильтр. Если же сначала сохранять в память отсчеты на всю длину фильтров, скажем N, то потребуется 2 умножения вместо 4 на каждый отсчет (за счет симметрии коэффициентов фильтров), но больше накладных расходов на организацию вычислений - сдается мне, экономии не получится, наоборот, зависит от процессора. При равной длине фильтров N, их коэффициенты займут в ROM 2N ячеек.

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

Таким образом вычисляется THD+N. Если нет шума, то оно будет равно THD.