subj
вообщем просто надоумьте мож я че не так понимаю.
сделал FFT звукового сигнала, построил спектр,
сделал отрисовку волны. все ок.
а как посчитать THD? из FFT это делается? или как то иначе.
сложить квадраты высших гармоник, поделить на квадрат основной гармоники, как то так вроде...
а реально как?

Сделать это можно согласно определению. Расчитываете мощность основной и высших гармоник, суммируете и делите. Можно это делать в частотной области через ДПФ, можно во временной пропуская сигнал через набор фильтров. Как вам больше нравится 

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

Надо вначале определится что это за сигнал. Основная гармоника может быть один отсчёт FFT а может несколько, зависит от того что подали на вход. Вы путаете спектральные отсчёты FFT и гармоники сигнала, это разные вещи. Для измерения THD удобно через систему пропусить гармонический сигнал и на нём проводить измерения. THD нельзя измерять на случайном сигнале, надо знать что подаётся на вход системы.

на входе синусоидный сигнал.
например 36Гц . дальше что?
это я читал http://en.wikipedia.org/wiki/Total_harmonic_distortion
чето не совсем понимаю как это запрограммировать...
получаю оцифрованные данные со зв. карты, делаю FFT, рисую спектр, рисую волну. а как из этих данных посчитать THD не понимаю чето...

Это сообщение отредактировал(а) pethead - 29.3.2011, 10:04

Если на входе 36 Гц то надо из сигнала выделить основную гармонику (36 Гц), вторую (72 Гц), третью (108 Гц) и так далее. Сделать это можно с помощью цифровых фильтров (это пожалуй будет самый простой вариант). Надо рассчитать фильтр нижних частот(ФНЧ) с частотой среза ~15...20 Гц и достаточным подавлением на отстройке 36 Гц (от этого зависит точность и диапазон измерения THD). Далее берёте свой сигнал и умножаете на комплексную экспоненту частоты 36 Гц, после умножения фильтруете сигнал, и на выходе получаем отсчёты амплитуды основной гармоники. Далее делаем те же операции но умножаем на частоты 72, 108 и т.д. Гц. Далее считаете мощности гармоник, и по формуле находим THD. Коэффициенты фильтра можно рассчитать в матлабе или аналогичной программе. Как получить отсчёты комплексной экспоненты и фильтровать сигнал надо писать ?

если не трудно напишите, но у меня цель все это сделать программно, без матлабов.
и еще: нельзя ли к расчету применить уже готовое преобразование фурье?
может есть где то готовый\полуготовый код по расчету THD? на любом языке, главное чтобы работал с оцифрованным сигналом. http://alglib.sources.ru/ в этом случае не может помочь?

Если хотите считать через спектр сигнала, то тогда так:

•  Рассчитываем размерность БПФ. Чтобы можно было использовать окно с хорошим подавлением боковых лепестков длину берём такой чтобы между гармониками попадало достаточно спектральных отсчётов БПФ. Предполагая что у вас частота дискретизации 44100 Гц, 2^16 на глаз будет достаточно.  
•  Делим сигнал на блоки длиной 2^16
•  На каждый блок накладываем окно (я попробовал с окном кайзера beta=14)
•  Считаем для каждого блока БПФ
•  Из полученных комплексных БПФ находим амплитудные спектры
•  Усредняем по нескольким блокам

Теперь на спектре нужно найти амплитуды гармоник частоты 36 Гц. Основной тон будет находится в районе 53 спектрального отсчёта 2-ая гармоника возле 106 и т.д. Но так как ваша частота не попадает точно в спектральные отсчёты то её максимум находится где то между ними. Максимум можно искать так:

•  Найти максимальный отсчёт на спектре
•  взять 2 соседние точки и построить по ним параболу
•  найти максимум параболы

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


Если хотите без расчёта фильтра то считайте через БПФ, как я написал.

Можно, как это сделать я написал выше. Только мне не ясно будет при этом выигрыш по скорости или проигрыш (считать неохота  ). Те операции что я описал в предыдущем посту через БПФ будут проделываться 2^16 раз вместо нескольких 10-ков раз. Смотря насколько сильные нелинейные искажения вашего тракта, и сколько придётся учитывать гармоник вашего сигнала. Если вам проще через БПФ считайте так. Принципиальной разницы нет.

Насчёт есть ли код это я не знаю, не интересовался. С alglib не работал.

Это сообщение отредактировал(а) vedun - 30.3.2011, 15:22

vedun, 
А если не через фильтры, а напрямую через гармонический анализ? Лучше же должно быть.





Да и вообще зачем фильтры и окна?


2^14 как-то понадежнее.
// Точность вычислений. При каждой операции мы теряем часть данных.
// Это часть маленькая, но так как элементов много она накапливается.
// Оценка дана в http://algolist.manual.ru/maths/long.zip
// Который ссылается на источник http://www.daemonology.net/tricl/fftconv_c.ps



Как по мне посчитать переходы через 0. Что позволит определить основную частоту более точно.  Заодно засечь начало первого и конец последнего периодов.

PS. Поправите меня, если я неправ.

Pavia насколько я разбираюсь в вопросе (гуру себя не считаю), то по сути это одно и тоже. Ведь алгоритм ДПФ это и есть последовательный перенос спектра с частоты Fd/N, 2*Fd/N, 3*Fd/N.... в 0, фильтрация ФНЧ (характеристика фильтра канала ДПФ вроде похожа на sin(x)/x но не она, какой там фильтр получается по памяти не скажу, но уровень АЧХ спадает достаточно медленно) и децимация в N раз. Зато если сразу работать по нужным частотам то не надо интерполировать спектр для поиска максимума.


Если не применять оконную функцию то надо брать большую размерность ДПФ чтобы избежать влияния гармоник друг на друга из за растекания спектра. Если к примеру уровень 2-й гармоники -80 дБ относительно 1-й то надо чтобы при измерении её уровня 1-ая была хотя бы на 10..15 дБ ниже 2-й, то есть подавление фильтра при отстройке 36 Гц должно быть 95 дБ. Я это на глазок по спектрам подобрал. Ещё непонятно насколько высок уровень нелинейных искажений у pethead и с какой точностью ему надо мерять THD, поэтому что либо советовать по окнам и размерности ДПФ сложно.

Я смотрел по спектрограммам на глаз, может хватит и 2^14.

Спасибо за ссылки, пока посмотрел по диагонали, как будет время изучу подробнее. Что могу сказать насчёт точности, построил спектр синусоиды с длиной ДПФ 2^24 степени (все вычисления с двойной точностью, тип double в С). Я не заметил чтобы ошибки сильно накопились (смотри рис.). Стоит одна честная палка 0 дБ остальные на уровне не более -320 дБ. По моему этого достаточно. Если в 2-х словах то как то так. 

vedun, Сам себя гуру не считаю. Вот и провоцирую народ, на дискуссию.

А вот тут можно подробнее? Где про это прочитать?

И это правильно, сам такой   

Например в книге Сергиенко "Цифровая обработка сигналов" стр. 270. 
Ну и классика жанра Рабинер, Голд. Сейчас вот туда заглянул (начиная со стр. 421) и обнаружил для себя новые вещи. Вот ещё статья о применении окон, на dsplib.ru вообще неплохая подборка статей по ДПФ и фильтрации.

спасибо. будем пробовать