анализ АЧХ звукового файла, или как это делается в демо.

Всем привет!
Собсно сабж. Не так давно делал демку со звуковым сопровождением. Нормально сделать программный фильтр по определённым частотным диапазонам так и не удалось. Хотелось бы узнать как это делается допустим в WinAmp и тому подобным программам, а делается там довольно качественно. И если есть теоретическая литература по поэтому вопросу просьба поделиться, (только не очень тяжёлая, а то у меня трафик очень маленький ).

 

Методов очень много, гугл в помощь. Ищи по преобразованию Фурье, КИХ и БИХ дискретным фильрам. Смотря насколько "качественным" нужно иметь результат. Есть много готовых примеров по Фурье, неплохая и быстрая либа есть тут http://fftw.org/

 

Bohdan200
Big thanks.
Очень высокое качество не требуется. В гугле пока ничего путного не нашёл . Но хоть теперь знаю куда копать .

 

Booster

1. Делаем FFT
2. Корректируем составляющие в нужном диапазоне.
3. Конвертируем обратно в PCM через IFFT.

Вот и всё. П1 и П3 самостоятельно реализовывать не нужно - есть уже готовые либы. На сайте Intel можно скачать либу с FFT и IFFT оптимизированную под Intel'овские процессоры. П2 - это набор простых арифметических действий. В общем, эквалайзер - штука очень простая.

 

Quantum
Штука она действительно довольно простая, но есть и свои "подводные камни"... Надо учитывать, что FFT предполагает, что входной сигнал периодический, на практике это не так и на выходном сигнале постле преобразования FFT->Умножаем на АЧХ->IFFT вблизи "крайних" (первого и последнего) значений появляются искажения в виде всплесков и т.п. Устраняется проделыванием операции с некоторым "нахлестом" и усреднения значений на этом "нахлесте" от двух состедних FFT

 

Bohdan200
Насколько я знаю, для этого используют взвешанную оконную функцию.
К примеру Хэмминга
(0.54-0.46*cos(2*pi*i/N))
Или это не так?

 

Pavia
Я брал ф-ю Гаусса

 

Quantum и Bohdan200
Можно поподробнее?
Пока пользовался fftw, на которую ссылался Bohdan200.
Насколько я понял, FFT преобразует график время/амплитуда, к графику частота/амплитуда.
На входе у нас N отсчётов амплитуд, а на выходе N/2 комплексных чисел. Первая составляющая этого числа - это амплитуда соответствующей гармоники, а вторая это фазовый сдвиг от начала(поправьте если не прав). Как я понимаю нужно вычислить текущую амплитуду соответствующей гармоники. Подозреваю что для этого нужно сделать что-то вроде sin(current_freq*t)+фаза - ?;
Пока пробовал брать небольшие временные интервалы, 64-1024 отсчётов и просто отображал амплитуду, то есть первую составляющую комплексного числа, но выходил такой же бардак, что и без использования библы fftw.

Quantum
Пожайлуста объясните поподробнее пункты 2 и 3 и что такое IFFT?. А также какое FFT нужно использовать, так как если я правильно понял есть разные подходы.

Bohdan200

Да это я увидел - первая гармоника и последня сильно искажены. Первая упирается в потолок, а последняя на нуле.
И еще, как я читал, в методе фурье амплитуды от больших частот к меньшим уменьшаются, и это на графике я увидел, но ведь в том же WinAmp, такого падения не наблюдается?
И еще на сайте библы fftw, примера для анализа АЧХ не нашёл (:, так что плиз, объясните более популярно, мне нужно сделать приемлимый спектроанализ в реальном времени, как это делается в WinAmp или Amarok. Но вообще хотелось бы сделать всё это самостоятельно, без либ, но для начала неплохо было бы сделать и с либами.

 

Booster
Использовать комплексное преоюразование имеет смысл только если интересует ФАЗА частотных составляющих. При этом АМПЛИТУДУ каждой из них можно подщитать по формуле
A=SQRT(R*R + I*I)
где R - действительная, а I - мнимая части комплексной амплитуды.
Если же фазовая компонента не интересует, можно использовать ДЕЙСТВИТЕЛЬНОЕ преобразование вместо полного комплексного. Выиграш по скорости - около 2-х раз. Действительное преобразование получает на входе и возвращает на выходе масив действительных отчетов. Подробнее можно почитать ТУТ, там очень неплохая статья с исходниками.
IFFT (Inverse Fast Fourier Transformation) - обратное преобразование Фурье.

Ничего никуда просто так не уменьшается, все зависит только от входного сигнала.

Чтото не так считаеш... Или входной сигнал имеет большую постоянную составляющую.

 

Bohdan200
Спасибо огромное.
Теперь должен с этим разобраться.
Еще решил немного помучать: A=SQRT(R*R + I*I) - это какая амплитуда, максимальная?. И что тогда есть R и I? Я думал что I это фаза гармоники, а R это коэффициэнт на который нужно умножить sin(f), поэтому думал, что для вычисления амплитуды нужно сделать A*sin(freq+I). Плиз разъясните чем же являются эти коэффициенты.

 

"Booster:
Использовать комплексное преобразование имеет смысл только если интересует ФАЗА частотных составляющих. "

Использовать комплексное преобразование имеет смысл почти всегда!
Если на входе действительный вектор, то в мнимую часть входных данных можно записать следующий "кадр" звуковых данных (или 2 канал, если стерео). Потом сделать комплексное БПФ - потом разделить частотные составляющие и получить спектр каждой функции. Таким образом можно одновременно вычислять 2 преобразования Фурье.
Ссылка: Э.Р. Канасевич "Анализ временных последовательностей в геофизике"
Если интересно могу подробнее рассказать...

 

2Santaev
Поподробнее было бы очень кстати.

2Bohdan200
На алголисте конечно статьи хорошие, но хотелось бы что-нить по теории в контексте обработки звука. Желательно чтобы формулами не сильно было перегружено.

Теперь пробовал по A=SQRT(R*R + I*I)
Беру примерно 1024-16386 отсчетов, делаю FFT, затем вычисляю среднее арифметическое по диапазонам частот, и отображаю скажем 16 полос. И график получается именно такой, что более высокие частоты имеют большие амплитуды. Может это потому что не использую взвешанную оконную функцию, или нужно использовать нахлёст?

 

Santaev

Я писАл, там по ссылке есть такой алгоритм. Вместо одного комплексного получаем два действительных. Выиграш около 2-х раз по скорости. Челу нужны просто "прыгающие столбики" как я понял...

Booster
Учти еже один момент. Чтобы все "красиво выглядело ", нужно зделать логарифмическую шкалу при отображении спектра. Слишком много отчотов лучше не бери, думаю 256 будет вполне достаточно. А оконная ф-я с нахлестом нужна только для "серьезных" применений FFT (DSP, эквалайзер...)

 

Bohdan200
думаю 256 будет вполне достаточно
Но тогда, слишком высока реакция, столбики слишком часто прыгают(мерцают) и это не красиво смотрится, да и проблему с разными высотами IMHO это не решает. А вот с логарифмической шкалой это интересно.

 

Знач так, алгоритм чисто для "столбиков" и чисто для "красоты"
* Делаем таймер (50мС к примеру);
* По таймеру делаем FFT того участка аудиобуфера, который в данный момент играет длинной 256 (128, 64, даже 32 для винампа хватит) отчета;
* Логарифмируем результаты (для "красоты" очень хорошо подходит SQRT);
* В проге должен быть некий буфер для усреднения результатов текущего FFT с результатами, скажем, 10 предыдущих. Усредняем;
* Рисуем результат.

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

 

Да исправят меня отцы.

Пример использования FFT:

Код (Text):

1. Cmplx getCmplx(double ampl, double phi){
2.     Cmplx res;
3.     res.Re = ampl * cos(phi);
4.     res.Im = ampl * sin(phi);
5.     return res;
6. }
7.  
8.  
9. double cmplxabs(Cmplx c){
10.     return sqrt(c.Re*c.Re + c.Im*c.Im);
11. }
12.  
13. #define SAMPLERATE 1024
14.  
15.  
16. int main(){
17.     Cmplx data[SAMPLERATE];
18.     double phi, ampl;
19.     for(int i = 0; i < SAMPLERATE; i++){
20.         phi = (double)i / SAMPLERATE * 2 * M_PI;
21.         // 3 составляющие: 16Гц, инт. 1, 8Гц, инт. 2, 128Гц, инт. 0.4
22.         ampl = sin(phi * 16 - M_PI / 3.0) + sin(phi * 8) * 2 + sin(phi * 500) * 0.4;
23.         data[i] = getCmplx(ampl, phi);
24.     }
25.     FFTC(data, SAMPLERATE);
26.  
27.     for(int i = SAMPLERATE / 2; i < SAMPLERATE; i++){
28.         if(cmplxabs(data[i]) > 1e-6){
29.             printf("%4d Hz: ampl = %lf; phase = %lf\n",
30.                 SAMPLERATE - i + 1,
31.                 sqrt(data[i].Im * data[i].Im + data[i].Re * data[i].Re) * 2 / SAMPLERATE,
32.                 atan(data[i].Im / data[i].Re));
33.         }
34.     }
35.  
36.     return 0;
37. }

 

32 выборки - маловато будет! Расстояние между частотными составляющими после БПФ будет равно Fd/N , где Fd -частота дискретизации , N- число отсчетов. При Fd=32 кГц между частотными составляющими будет 500 Гц.
Так что в области низких частот отсчетов очень мало. Звуковой сигнал лучше рассматривать по октавам : например 10-20, 20-40, 40-80 .... 10240-20480 Гц, соответственно и величину N надо выбирать для каждого диапазона частот.
Т.е для "красоты" логарифмировать надо и по амплитуде, и по частоте!
По амплитуде - чтобы "визуализировать" широкий динамический диапазон звукового сигнала (~168 дБ), но если у тебя 8-битные данные (~48 дБ), то SQRT тоже хорошо!

 

Я извиняюсь - при N=32 и Fd=32кГц между частотными составляющими -1кГц!

 

2Santaev
Big thanks too.
Действительно выбирать различные N и логарифмировать по частоте - отличная идея. А то у меня при 16 равномерых полос, в основном только средние и были видны
Юзаю 16-бит данные, как я понимаю 96 дб.
Интересно а логарифмировать лучше по какому основанию?
Частоту как я понимаю по 2. А амплитуду пробовал по 10, вроде не плохо, но может слудует по другому?
И ещё, первую частоту я вообще не показываю так как её амплитуда вообще почти максимальна. Почему это может быть?

 

Это будет, грубо говоря, зависеть от количества пикселей которые ты отведешь на визуализацию. (или 96 дБ "спрессовать" в 10 пикселей или в 1024 - будет разница). Это не проблема подобрать.

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