Преобразование фурье

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

2СанитарЖеня

непонял???????

Вот пример реализации на Паскале. Здесь все прозрачно.

procedure fft(var Ar:array of double; var Ai:array of double);
// На входе - действительная и мнимая части, на выходе они же, но для образа
var
Le,Le1,jt,it,ip,lt,kt:Integer;
Tr,Ti,Ur,Ui,Wr,Wi:Double;
// Предполагается, что точек 256, отсюда константы 255, 128 и 8.
// Изменения и обобщения очевидны
begin
    for it:=0 to 255 do
    begin
    kt:=128;
    jt:=0;
    Le:=1;
    Le1:=it;
    for lt:=1 to 8 do
    begin
     if kt<=Le1 then begin
     jt:=jt+Le;
     Le1:=Le1-kt;
     end;
     Le:=Le*2;
     kt:=kt div 2;
    end;
// Это предобработка - перестановка элементов в порядке инвертированных битов №
    if it<jt then begin
       Tr:=Ar[jt];Ti:=Ai[jt];
       Ar[jt]:=Ar[it];Ai[jt]:=Ai[it];
       Ar[it]:=Tr;Ai[it]:=Ti;
       end;
    end;

// Основной цикл
    for lt:=1 to 8 do
    begin
     Le:=1;
     for jt:=1 to lt do Le:=Le*2;
     Le1:=Le div 2;
     Ur:=1.0;
     Ui:=0.0;
// Значение синуса и косинуса вычисляется по формуле С. и К. суммы углов
     Wr:=cos(Pi/Le1);
     Wi:=sin(Pi/Le1);
// Проход по элементам, с постоянно уменьшающимся шагом
     for jt:=0 to Le1-1 do
     begin
      it:=jt;
      while it<256-Le1 do
      begin
       ip:=it+Le1;
// Бабочка. Оптимизирована по скорости
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;
       it:=it+Le;
      end;
// Вычисление новых значений синуса и косинуса
     Tr:=Ur*Wr-Ui*Wi;
     Ti:=Ur*Wi+Ui*Wr;
     Ur:=Tr;
     Ui:=Ti;
     end;
    end;
end;

2СанитарЖеня

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

Эээ... А понять ответ лень не помешает?
Поворачивающие множители - связаны с известным свойством Фурье, а именно, что при сдвиге (а объединяемые отрезки сдвинуты относительно друг друга) Фурье-образ претерпевает поворот (по фазе сдвигается, попросту).
А откуда новые частоты - оттуда, где и были. Но, для короткого отрезка, они не могли быть выделены явно, а искажали образы отрезков. По разному искажали. Сумма образов отрезков (поправленных корректирующими множителями) дает уточненную величину тех частот, которые были явно видимы в образах, а разность - дает это искажение, по которому и восстанавливается образ в промежуточных точках.

2СанитарЖеня

сколько ни силился, всё равно не понял. твоё объяснение похоже на древний алхимический текст. ты бы чтоль формулу написал, как выходят коэффициенты при нечётных гармониках?

Все к Рабинеру и Голду (или, кто не желает - к Рабинеру и Гоулду). dsp-book.narod.ru
Там подробно изложено. Можно и Гольденберга, Матюшкина и Поляка.
Вкратце - вычитаем две оценки четных компонент, и их разность дает нечетную компоненту. Ну и те же поворачивающие множители.

2СанитарЖеня

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

Ну вот же они, формулы. Вот она, бабочка. Что именно неясно?        
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;

Цитата СанитарЖеня, 31.03.03, 11:20:58

Ну вот же они, формулы. Вот она, бабочка. Что именно неясно?        
       Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
       Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
       Ar[ip]:=Ar[it]-Tr;
       Ai[ip]:=Ai[it]-Ti;
       Ar[it]:=Ar[it]+Tr;
       Ai[it]:=Ai[it]+Ti;

я ща заплачу (заплАчу :'(). чтотакое tr, ar, ur, далее везде???????????????????

Господи... Ну как можно объяснить проще? Ну нет в ДСП объяснений на уровне "Ты мне пальцем ткни!" - считать надо, и думать...
Кристально ясный же алгоритм...
Ну почитайте же Рабинера и Голда (АКА Гоулд).
dsp-book.narod.ru
Или Гольденберга, Матюшкина и Поляка
dspbook.km.ru.

Что есть Ar и Ai - в комментарии вначале программы написано, Ur и Ui - действительная и мнимая часть поворачивающего множителя, Tr и Ti - промежуточные переменные, необходимость их вызвана тем, что операции выполняются последовательно.

Очень нужна помощь в организации Фурье-фильтра в Matlab.
Смысл задачи:имеется 2 массива данных, один - частоты, второй - поглощение по ним (фурье-спектроскопия веществ). Необходимо, выделив одну частоту получить, на выходе поглощения по этой частоте.
Заранее благодарен, это очень сторочно!
P.S. Сам я в Матлабе ламер, если есть возможность, то опишите подробно.

to Ciceron
Нужно делать цифровой фильтр (полосовой), по-моему в MatLab есть примеры.

to All
Есть несколько приемов получения FFT (sic!) для призвольного числа точек (<>2^N).
Самый простый -- дополнение 0-ми до 2^N (перерасход памяти и снижение скорости). Посложнее -- разложение на простые множители для, которые и вычисляются аналитически, в случае =2^N этим простым множителем является 2 -->"бабочка", аналогично можно расписать преобразование для 3, 5, 7 и  т.д. Самый   худший случай когда число точек простое и необходимо выполнять DFT.

Чуваки, дайте лучше сцылку на оптимизированное 16бит в 16ит преобразование без использования FPU

Разобрался

Сообщение отредактировано: greek - 03.01.06, 19:00