%fft of x
%x is the input vector
%N is the number of the points
%at 024009 2009-4-3
function y = fft_custom(x,N)
%逆序操作
temp_reverse=[0:N-1];
LH=N/2;
J=LH;
N1=log2(N);
for I=temp_reverse(2):temp_reverse(N-1)
if I<J
T=temp_reverse(I+1);
temp_reverse(I+1)=temp_reverse(J+1);
temp_reverse(J+1)=T;
end
K=LH;
while J>=K
J=J-K;
K=K/2;
end
J=J+K;
end
temp_reverse=temp_reverse+1;
for i=1:N
y(i)=x(temp_reverse(i));
end
%蝶形计算
M=log2(N);
Wn=exp^(-j*2*pi/N);
for L=1:M %表示蝶形运算的级数
B=2^(L-1); %每个蝶形输入数据相距的点数,和不同蝶形数目
for J=0:B-1; %每级参与计算的蝶形数
P=2^(M-L)*J; %当前蝶形中旋转因子的指数
for k=J:2^L:N-1;
T=y(k+1)+y(k+1+B)*Wn^P;
y(k+1+B)=y(k+1)-y(k+1+B)*Wn^P;
y(k+1)=T;
end
end
end
评论