参考文献《数字信号处理. 理论、算法与实现》 第二版实现DIT基2FFT,主要就是实现下图中的一个蝶形运算。
下面是用matlab编写的实现函数。
function [Xk]=DIT_FFT_2_MOD(xn,N)
t=1:N;
WWr=exp(-1i*2*pi/N*[0:N/2-1].');%旋转因子
%蝶形运算开始
M=log2(N);%“级”的数量
% 码位倒置
Xk=permute(reshape(xn,2*ones(1,M)),M:-1:1);
Xk=Xk(:);
N2=N/2;
Num_of_Group=N2;%每一级中组的个数初始值
Interval_of_Group=1;%每一级中组与组之间的间距
for m=0:M-1 %“级”循环开始
Wr=WWr(1:Num_of_Group:end);
gMatrix2=reshape(t,Interval_of_Group,2,Num_of_Group);
gMatrix21=reshape(gMatrix2(:,1,:),N2,1);
gMatrix22=reshape(gMatrix2(:,2,:),N2,1);
if(m==0)
XKtemp=Xk(gMatrix22);%第0级,蝶形运算式
elseif(m==1)
XKtemp=Xk(gMatrix22);
XKtemp(2:2:end)=XKtemp(2:2:end)*Wr(2);%第1级,蝶形运算式
else
ss=repmat(Wr,Num_of_Group,1);
XKtemp=Xk(gMatrix22).*ss;%第m级,第g组的蝶形运算式1
end
Xk(gMatrix22)=Xk(gMatrix21)-XKtemp;%第m级,蝶形运算式
Xk(gMatrix21)=Xk(gMatrix21)+XKtemp;
Interval_of_Group=Interval_of_Group*2; %递推
Num_of_Group=Num_of_Group/2; %递推
end
end
下面调用函数并测试16点FFT的运算。
clc
N=16; % 数据长度为2的整数幂次
x=(0:N-1).'; %变换的数据,列向量格式
Xf=DIT_FFT_2_MOD(x,N); %FFT变换,长度N在反变换除
Xk=fft(x);
[Xf.';Xk.']
下面是一次运行的结果
ans =
1.0e+02 *
Columns 1 through 7
1.2000 + 0.0000i -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i -0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i
1.2000 + 0.0000i -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i -0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i
Columns 8 through 14
-0.0800 + 0.0159i -0.0800 + 0.0000i -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i -0.0800 - 0.0800i -0.0800 - 0.1197i
-0.0800 + 0.0159i -0.0800 + 0.0000i -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i -0.0800 - 0.0800i -0.0800 - 0.1197i
Columns 15 through 16
-0.0800 - 0.1931i -0.0800 - 0.4022i
-0.0800 - 0.1931i -0.0800 - 0.4022i