- 空间
function X = fft_dit2(x,f)
% 时域抽取基2-FFT(基2-DIT-FFT)
% 刘顺兰,吴杰,数字信号处理,Page 113.
% ----------------------------------------
% function X = fft_dit2(x)
% x = 输入数组(长度为2的幂)
% f = 1:Inverse, -1:Forward
% X = 输出数组
%
% isReal = isreal(x);
nx = length(x);
M = 0;
N = 1;
while N < nx
M = M+1;
N = 2*N;
end
X = [x(:);zeros(N-nx,1)];
% =================================
% 倒位序
% X = x;
LH = N/2;
J = LH;
N1 = N-2;
for I = 1:N1
if I < J
T = X(I+1);
X(I+1) = X(J+1);
X(J+1) = T;
end
K = LH;
while J >= K
J = J-K;
K = K/2;
end
J = J+K;
end
% ==============================
C = -j*2*pi*(-f)/N;
for m = 1:M % 共M级,最外层循环M次
B = 2^(m-1); % B = 第m级的蝶形两点间距
for J = 0:B-1 %2^(m-1)-1 % 第m级有2^(m-1)个不同的旋转因子W^J,J = 0,1,...,2^(m-1)-1
P = 2^(M-m)*J;
Wnp = exp(C*P); % 旋转因子
% gn = 2^(M-m) % 第m级,组数gn
for k = J:2^m:N-1 % 每个旋转因子对应2^(M-m)组蝶形计算
tmp1 = X(k+1); % 蝶形计算
tmp2 = X(k+B+1)*Wnp;
X(k+1) = tmp1+tmp2;
X(k+B+1) = tmp1-tmp2; % 同址计算
end
end
end
% Inverse Transform
if f == 1
X = X/N;
end
收藏于 2013-08-28
来自于百度空间