Отправлено Vadim Kudryavtsev 17 января 2005 г. 15:43 В ответ на: Ответ: По книжке Блейхута "Быстрые алгоритмы цифровой обработки сигналов"(+) отправлено Vadim Kudryavtsev 15 января 2005 г. 09:00

Выдалось полтора часа свободных - сотворил для собственного интереса быстрый алгоритм для длины 50. Интересно, что в результате используются только преобразования длины 2 и 4.
Итого:
1. Преобразование длины 50 разбивается на 25 преобразований длины 2 и 2 преобразования длины 25.
2. Каждое из преобразований длины 25 разбивается на 5 преобразований длины 5 и еще 5 преобразований длины 5 (те реализуется десятью преобразованиями длины 5)
3. Каждое из преобразований длины 5 вычисляется с помощью трех преобразований длины 4.

Итого - 25 преобразований длины 2 (просто бабочек), и 2*10*3=60 преобразований длины 4 (в каждом 4 бабочки и умножение на j)
На шаге 2 требуется 2*16=32 умножения на не тривиальные комплексные числа.
На шаге 3 - 2*10*4=80 умножений на не тривиальные комплексные числа.

Правильно Блейхут писал - практически для любой длины можно найти эффективный алгоритм.

Остается при реализации продумать, как грамотно поступать с пред и пост перестановками данных

Код на МатЛабе

function y=fft50(x);

% for debug
% x=randn(1,50);

% mapping "x" into square matrix
s=zeros(2,25);
for i=0:50-1, s(rem(i,2)+1,rem(i,25)+1)=x(i+1); end

% apply fft to each column
s1=fft(s); % 50 fft's of length 2 (butterflies)

% apply fft to each row
%ss=conj(fft(conj(s1'))');
for i=1:2, ss(i,:)=fft25(s1(i,:)); end % 2 fft's of length 25

% map to output vector
y=zeros(1,25);
for i=0:50-1, y(rem(rem(i,2)*25+rem(i,25)*2,50)+1)=ss(rem(i,2)+1,rem(i*1,25)+1); end

function y=fft25(x);

% for debug
% x=randn(1,25);

% mapping "x" into square matrix
s=zeros(5);
for i=1:5, s(i,:)=x((i-1)*5+1:(i-1)*5+5); end

% apply fft to each column
% s1=fft5(s);
s1=zeros(5);
for i=1:5, s1(:,i)=conj(fft5(conj(s(:,i)'))'); end % 5 fft's of length 5 (transpose simultaneously conjugates vector)

% multiply by twidle factors
i=sqrt(-1);
s2=s1.*exp(-i*2*pi/25*[0:4]'*[0:4]);

% apply fft to each row
%ss=conj(fft5(conj(s2'))');
ss=zeros(5);
for i=1:5, ss(i,:)=fft5(s2(i,:)); end % 5 fft's of length 5

% map to output vector
y=zeros(1,25);
for i=1:5, y((i-1)*5+1:(i-1)*5+5)=ss(:,i); end

function y=fft5(x);

% for debug
% x=randn(1,5);
% fft(x)

i=sqrt(-1);
w=exp(-i*2*pi/5*[0:4]);

input_perm=[1 2 3 5 4]; % input permutation law
output_perm=[1 3 2 4 5]; % output permutation law

x=x(input_perm);
w=w(output_perm);

y=zeros(1,5);

y(2:5)=ifft(fft(w(2:5)).*fft((x(2:5)))); % 3 fft's has length 4

y(1)=sum(x);
y(2:5)=y(2:5)+x(1);

y(output_perm)=y;