4. Алгоритм быстрого преобразования Фурье

Запишем N-точечное дискретное преобразование Фурье для входной последовательности X_m.

. (1)

- 123 -

Разделим входную последовательность X_m отсчетных значений на четные A_m и нечетные B_m отсчеты

. (2)

Обозначим через и k-е компоненты преобразованных по ДПФ последовательностей A_i и B_i.. Тогда (2) можно записать более компактно

. (3)

Из (2-3) видно, что одно N-точечное ДПФ можно заменить двумя N/2-точечными ДПФ. Такая замена позволяет существенно уменьшить количество операций, затрачиваемых на вычисление всех компонент ДПФ. Рис.10 иллюстрирует процедуру вычисления ДПФ, получившую название «прореживание по времени». В этой процедуре выполняется замена одного 8-точечного ДПФ двумя 4-точечными ДПФ. Процедура замены результатов действия двух N/2-точечных на N-точечное ДПФ согласно выражению (3) получила в литературе название «бабочка». На рис.10 эта операция обозначена в виде точек с заливкой в пересечении соответствующих линий.

Рис. 10. Алгоритм вычисления ДПФ с прореживанием по времени. Реализация 8-точечного ДПФ с помощью двух 4-точечных ДПФ.

Оценим выигрыш по количеству операций, затрачиваемых на вычисление ДПФ различными способами. Для N-точечного ДПФ необходимо выполнить N² операций умножения отсчетных значений на значения

- 124 -

соответствующей гармонической функции. Для N/2-точечного ДПФ необходимо выполнить 2*(N/2)² операций умножения отсчетных значений на значения соответствующих гармонических функций и N/2 операций умножения преобразованных значений на экспоненту. Общее сокращение количества операций составит раза,

т.е. примерно в 2 раза. Если N является целой степенью числа 2, то такой процесс замены количества точек в ДПФ можно продолжить и далее до 2-точечного ДПФ включительно. В 2-точечном ДПФ не требуется умножения отсчетных значений на гармоническую функцию, поскольку последняя принимает лишь два возможных значения +1 и –1.За счет этого в 2-точечном ДПФ необходимо нахождение лишь суммы и разности двух отсчетных значений. Соответствующий алгоритм, использующий такую последовательную замену на ДПФ с меньшим количеством точек до 2-точечного ДПФ, получил название быстрого преобразования Фурье (БПФ). Общий выигрыш по количеству операций при использовании БПФ составляет

раз.

Например, если N=2048, то объем вычислений при использовании БПФ, по сравнению с ДПФ сократиться в

раз.

Рассмотрим текст процедуры БПФ, которая может быть использована при работе в среде Турбо-Паскаля или Delphi.

1) Процедура БПФ состоит из двух частей: собственно преобразования Фурье и циклов перестановок преобразованных значений x[i], y[i] в порядке возрастания частоты.

2) Входной процесс и результат БПФ в комплексной форме записываются в одни и те же массивы x[i], y[i]. БПФ выполняется в трех вложенных циклах. Двухточечное комплексное ДПФ осуществляется во внутреннем цикле.

3) Объем выборки N , обозначенный через Num, обязательно должен быть равен 2^MM1, где MM1 - целое число.

4) Расчет значений синуса и косинуса, необходимых для выполнения вычислений по формуле (3), осуществляется по рекуррентным формулам, аналогичным рассмотренным ранее. При этом s=sin(2/N), c=cos(2/N). Текущее значение косинуса заносится в U1, синуса - в U2.

5) Описанная процедура БПФ позволяет вычислять как прямое (inv=1), так и обратное (inv=-1) преобразование Фурье.

- 125 -

Текст процедуры быстрого преобразования Фурье

procedure FFT(Nn,inv:integer);

var

t1,t2,U1,U2,U3,C : real;

MM1,L0,L,L1,J,i1,N1,N2,k : integer;

Begin

MM1:=round(ln(Nn)/ln(2)); L0:=Num;

{------------------------ Блок БПФ -------------------------}

for L:=1 to MM1 do begin

L1:=L0 div 2 ;

U1:=1 ; U2:=0 ;

C:=cos(pi/L1) ;

S:=inv*sin(pi/L1) ;

for j:=1 to L1 do begin

{---------------- Начало внутреннего цикла -----------------}

i:=j; while i<=Num do begin i1:=i+L1;

T1:=x[i]-x[i1] ; T2:=y[i]-y[i1];

x[i]:=x[i]+x[i1] ; y[i]:=y[i]+y[i1];

x[i1]:=T1*U1-T2*U2 ; y[i1]:=T2*U1+T1*U2;

i:=i+L0 ; end;

u3:=U1*c-U2*s ; U2:=U2*c+U1*s ; U1:=U3 ; end;

L0:=L0 div 2 ; end;

{----------------- Перестановки -------------------------}

N2:=Num div 2; N1:=Num-1; j:=1;

for i:=1 to N1 do begin

if i<j then begin

T1:=x[j] ; T2:=y[j] ;

x[j]:=x[i] ; y[j]:=y[i] ;

x[i]:=T1 ; y[i]:=T2 ;

end; k:=N2; while k<j do begin j:=j-k ; k:=k div 2 ; end;

j:=j+k ; end ;

end;

• 126 –