Вот fft.h: Код (Text): // Быстрое преобразование Фурье #ifndef _FFT_ #define _FFT_ // Прямое быстрое преобразование Фурье. // Действительная и мнимая части отсчетов исходного сигнала // должны быть помещены в массивы a и b (одинаковой длины). // Длина массивов (n) должна быть целой степенью (r) двойки. // После выхода из функции эти же массивы содержат // вещественную и мнимую части спектров. void cfft( double *a, double *b, int r, int n ); // Обратное быстрое преобразование Фурье. // Действительная и мнимая части отсчетов спектра // должны быть помещены в массивы a и b (одинаковой длины). // Длина массивов (n) должна быть целой степенью (r) двойки. // После выхода из функции эти же массивы содержат // вещественную и мнимую части соответствующего спектру сигнала. void cifft( double *a, double *b, int r, int n ); #endif // _FFT_ Вот fft.cpp: Код (Text): #include "fft.h" #include <math.h> #define FFT_PI 3.14159265358979323846 static void binrevers( double *a, int n ) { int i, j, k; for( j = i = 0; i < n - 1; i++ ){ if( i < j ){ double t = a[j]; a[j] = a[i]; a[i] = t; } for( k = n / 2; k <= j; k /= 2 ) j -= k; j += k; } } static void base_fft( double *x, double *y, int r, int n ) { double ur, ui, wr, wi, tr, ti, ur2; int i, j, l, le1, le2, ip; binrevers(x, n); binrevers(y, n); for( le2 = l = 1; l <= r; l++ ){ le1 = le2; le2 *= 2; ur = 1.0; ui = 0.0; wr = cos(FFT_PI / le1); wi = -sin(FFT_PI / le1); for( j = 0; j < le1; j++ ){ for( i = j; i < n; i += le2 ){ ip = i + le1; tr = x[ip] * ur - y[ip] * ui; ti = x[ip] * ui + y[ip] * ur; x[ip] = x[i]; y[ip] = y[i]; x[i] += tr; y[i] += ti; x[ip] -= tr; y[ip] -= ti; } ur2 = ur * wr - ui * wi; ui = ur * wi + ui * wr; ur = ur2; } } } static void scale( double *a, int n ) { int i; for( i = 0; i < n; i++ ) a[i] /= n; } static void resort( double *a, int n ) { int i; for( i = 1; i < n / 2; i++ ){ double t = a[i]; a[i] = a[n - i]; a[n - i] = t; } } void cfft( double *a, double *b, int r, int n ) { scale(a, n); scale(b, n); base_fft(a, b, r, n); } void cifft( double *a, double *b, int r, int n ) { resort(a, n); resort(b, n); base_fft(a, b, r, n); } Это самый маленький рабочий сорец, который я нашел в инете. Никак не могу понять, почему восстановление функции по спектру происходит тем же извращенным способом, что и вычисление спектра? Что за спектр возвращает FFT?
_DEN_ Вопрос по теории или "извращенной" реализации ? Если по теории, то прямое и обратное преобразование отличаются постоянным множителем 1/n и изменением знака аргумента в sin\cos. В приведенном сорце это учитывается вызовами scale и resort. "Что за спектр возвращает FFT" - вопрос странный, если ты ищешь его реализацию и не знаешь что это такое. FFT возвращает дискретный спектр амплитуд a и b функций cos и sin на частотах 2*pi*i/n (i = 0..n-1). Если просуммировать эти функции, то получим тригонометрическую интерполяцию исходного сигнала на интервале 0..n-1 (и периодически повторяющуюся вне этого интервала). В заданных точках 0..n-1 эта сумма точно совпадает с отсчетами сигнала. Ну я думаю ты и сам представляешь. Что касается "извращенности", да шут ее знает. На первый взгляд ничего, а подробнее копать лень. PS: что мне лично не по душе, это когда раздельно вызывают sin и cos одного аргумента. Чего там компилятор наоптимизирует - неизвестно. Если math.h от борланда, то там точно должна быть процедурка SinCos, использующая инструкцию FPU FSinCos - в полтора раза быстрее чем по отдельности.
> Чего там компилятор наоптимизирует - неизвестно Нормальный компилятор все нормально соптимизирует даже без внешних либ, вот например шо делает интел. .text:000000AE loc_AE: ; CODE XREF: base_fft(double *,double *,int,int)+142j .text:000000AE fld ds:_2il0floatpacket_1 .text:000000B4 fld ds:_2il0floatpacket_3 .text:000000BA mov esi, edx .text:000000BC mov [esp+30h+var_30], esi .text:000000BF fild [esp+30h+var_30] .text:000000C2 fdivr st, st(3) .text:000000C4 add edx, edx .text:000000C6 fsincos
leo Приведи пример для данного сорца. На вход даем: a[2] = {4, 7}; b[2] = {0, 0}; Получаем спектр: a[2] = {5.5, -1.5}; b[2] = {0, 0}; Вот и коим хреном без извращений восстановить функцию? Тоесть хотелось бы увидеть просто сумму косинусов / синусов. Именно для этих двух чисел. И еще не совсем ясно, как ФФТ представляет себе дискретный кусок, четной функцией или нечетной?
_DEN_ Если ты под "извращением" понимаешь не конкретную реализацию, а метод БПФ вообще, то ты не прав. При расчете ДПФ в лоб, мы имеем большое число повторяющихся операций умножения отсчета сигнала на одно и тоже значение sin\cos т.к. a) это периодические функции, б) с кратными периодами, в) sin(x) = cos(x+pi/2). Метод БПФ уменьшает число таких операций грубо с n^2 до n*log2(n) и достигается это за счет "извращения": особой группировки отсчетов сигнала и вычисления спектра не последовательно отсчет за отсчетом, а хитро - на каждом шаге получаем вклад в два комплексных отсчета спектра. "хотелось бы увидеть просто сумму косинусов / синусов" Пусть отсчеты сигнала это S=x+j*y, а спектр F=a+j*b, тогда S(i)=Сумма(F[k]*exp(j*2*pi*i*k/n)), k = 0..n-1 Для действительной части получим: x(i)=Сумма(a[k]*cos(2*pi*i*k/n)-b[k]*sin(2*pi*i*k/n)), k = 0..n-1 для твоего примера: x[0]=5.5-0-1.5-0 = 4 x[1]=5.5-0-1.5*cos(2*pi*1*1/2)-0 = 5.5+1.5 = 7 "четной функцией или нечетной ?" Строго говоря, спектр ДПФ соответствует сигналу S(i) периодически продолженному в область i < 0 и i >= n. Поэтому четность\нечетность определяется видом S(i): если S(i) симметрична относительно n/2 - при продолжении в i < 0 получим четную функцию, если антисимметрична - нечетную, если вообще несимметрична, то ... Хотя вся эта четность\нечетность - по сути одна казуистика. PS: если тебе не нужно быстродействие или нужен просчет спектра не во всех точках, можешь плюнуть на этот БПФ и считать в лоб.
Dr.Golova "Нормальный компилятор все нормально соптимизирует" Я и раньше сомневался во всемогуществе "нормальных компиляторов", а nOp & volodya окончательно убедили меня в этом: "В этом языке НЕТ ограничений. Это всякие высокоуровневые языки следят за программером, как за дитем. В асме все по-взрослому. Суровые ребята пишут суровый код." (С) 2004, nOp PS: в каждой шутке есть доля истины
leo Нет, FFT в целом извращением я не считаю. Снижение с o(n^2) до o(log(n,2)*n) это рулез. На данном примере через простую сумму ряда все действительно работает. Но вот загрузил wav-ку с голосом, разложил по fft и свернул обратно через простое суммирование ряда - и нифига не получилось В наушниках какое-то монотонное журчание Что я в итоге хочу. У меня есть wav-ка. Это либо голос, либо какой-то семпл. Я хотел попробовать получить его спектр через sin-cos, а восстановить по какой-нибудь другой ортонормированной последовательности. Например sin^3(x), cos^3(x). Нормированной она конечно не будет, но ортогональной будет. Но в принципе - нормировать ортогональную последовательность не проблемма - лишь бы первообразные произведения двух произвольных членов выражались через элементарные функции. Ладно, хрен с ним Что вообще можно сделать с сигналом (читайте - "со звуком"), если мы знаем его дискретный спектр и можем с ним поизвращаться, но восстанавливать будем только по sin-cos? Ну, понятное дело, что первое что приходит на ум это частотные фильтры. Но это не очень интересно Что хорошего из спектра еще можно "выжать"?
>В наушниках какое-то монотонное журчание что бы не было журчания, нужно раскладывать не целиком wav-ик, а разделить его на достаточно малые "куски", и раскладывать, и собирать каждый по отдельности. >Что хорошего из спектра еще можно "выжать"? ну в основном это для фильтров и применяется, так же его можно сжать, просеивая через психо-аккустический, опять же, фильтр. (Впрочем это и есть MP3 сжатие)
Johnikum А какая хрен разница-то? Или там с ограничением мантиссы хрень происходит? Точности что-ли не хватает. Да мне не надо его сжимать - я хочу получить какие-нибудь эффекты, дешевые по размеру... Например, мне кажется (догадка), что если вместо cos(x) взять pow(abs(cos(x)),1+k/(t+1))*sign(cos(x)), то получится плавный переход от нормального звучания к железно-кислотному... Хотя это только догадка...
А какая хрен разница-то? Вообще говоря, перевод сигнала в частотное представление возможен только блоками (окнами). В нашем случае - FFT блоками. Мы делим исходный сигнал на блоки и говорим: в этом блоке имеются такие-то частоты. Вернее, так: его можно получить, сложив такие-то частоты. FFT раскладывает функцию не на её гармоники, а на свои гармоники. Допустим, в исходной функции была лишь частота точно 100 Гц, а мы с помощью FFT раскладываем функцию в частоты 96, 99, 102, 105 Гц - как вы думаете, что получится в результате? Ничего полезного. FFT по прежнему будет способна полностью точно восстановить исходную функцию, однако полученный нами спектральный состав будет бестолковым. взято отсюда: www.mtu-net.ru документ называется: m_fft1.htm
Ну уж совсем за дебила-то меня держать не надо Я этот cpp-шник за 5 минут на асм переложу Вот! Как же тогда из FFT-шного спектра получить нормальный?
"не на её гармоники, а на свои гармоники" "из FFT-шного спектра получить нормальный" Пример со 100 Гц это общая проблема расчета спектра численными методами. При любом методе мы должны задать дискрет расчета по частоте fmin и всегда в сигнале может найтись составляющая, которая попадет между этими дискретами (если возьмем 1 Гц, то "провалится" 100.5 Гц и т.д.). "Свои гармоники" это i*fmin. В этом смысле особенностью FFT является только то, что fmin = 1/T = 1/(n*t), где n = 2^r (t - интервал времени между отсчетами, T = n*t - длительность анализируемого блока). Ясно, что если 100 Гц попадает между дискретами i*fmin (т.е. не кратна частоте дискретизации сигнала), то расчетный спектр "размазывается" и тем больше, чем меньше n. А почему - грубо говоря потому, что не выполняется условие ортогональности, т.е. на малом интервале при малом числе отсчетов почти все суммы произведений cos 100 Гц и cos i*fmin оказываются не равными нулю. По сути выход из этой ситуации только в увеличении n. Но т.к. при увеличении n ухудшается динамика, то используют разнообразные весовые окна. (Можно также трактовать действие окна как раширение спектра исходного сигнала, т.е. вместо синусоиды 100 Гц мы получим после умножения на окно импульсный сигнал, спектр которого занимает несколько наших дискретов i*fmin, соответственно всплеск на этих частотах усиливается, а вклад далеких уменьшается).
_DEN_ Ну уж совсем за дебила-то меня держать не надо Извините, но даже не было в мыслях обижать. Вот! Как же тогда из FFT-шного спектра получить нормальный? Идеально точно никак.
FFT, оно же БПФ, оно же вариант дискретного преобразования Фурье с оптимизацией вычислений, поэтому спектр, полученный после преобразования FFT есть дискретный спектр. По поводу размытия спектра (to leo) в случае если частота сигнала совпадает с i*df, где df = частота_дискретизации/длина_выборки, то будет один максимум, если совпадает с (i+0.5)*df, то два одинаковых максимума. Ширина лепестка (ненулевых составляющих вокруг максимумов) зависит от того, целое ли число периодов сигнала входит в выборку. Если целое, то ширина лепестка минимальная, иначе она увеличивается. Для борьбы с этим явлением используют оконные функции. В двух словах все
_DEN_ спрашивал: Что вообще можно сделать с сигналом (читайте - "со звуком"), если мы знаем его дискретный спектр и можем с ним поизвращаться, но восстанавливать будем только по sin-cos? Ну, понятное дело, что первое что приходит на ум это частотные фильтры. Но это не очень интересно Что хорошего из спектра еще можно "выжать"? javascript:insert_text('[bb]\n', '') Очень БПФ эффективно при вычислении функции взаимной корреляции... На практике распознаватель речи можно построить простой (хотя бы на цифры)... Насчет примера численного у меня некоторые сомнения:в спектре b[2] = {0, 0}; ??? Для действительной функции на входе вектор а должен быть четной функцией, а b - нечетной...Причем в первой половине а ,b - положительные частоты - а во второй отрицательные. Странно также, что при обратном преобразовании получился не такой же звук. Функция при обратном преобразовании восстанавливается 100%. Может wav-ик неправильно раскрыл? (там же вроде целые без знака) или при преобразовании действительных в wav ошибка? Да, кстати, скоро юбилей - преобразованию Фурье (только не быстрому)- 200 лет !
