DIT和DIF实现快速傅里叶变换的FFT
数字信号处理课程老师要求做的一个小程序,在此分享~
%***************基2的DIF&DIT-FFT算法************************
%****************——徐**——*****************************
%
%*****************分组:第三组******************************
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%频率f1、f2输入
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YN='Y';
while(strcmp(YN'Y'))
f1=input('请输入频率f1的值(200<=f1<300):');
while((f1<200) || (f1>=300))
disp('输入的频率不满足"200<=f1<300"的条件,请重新输入:');
f1=input('f1:');
end
f2=input('请输入频率f2的值(300<=f2<400f1+f2>600):');
while((f2<300) || (f2>=400))
disp('输入的频率不满足"300<=f2<400"的条件,请重新输入:');
f2=input('f2:');
end
while (f1+f2<=600) %判断输入的f1、f2是否满足条件
disp('输入的频率不满足"f1+f2>600"的条件,请重新输入:');
f1=input('请输入频率f1的值(200<=f1<300):');
while((f1<200) || (f1>=300))
disp('输入的频率不满足"200<=f1<300"的条件,请重新输入:');
f1=input('f1:');
end
f2=input('请输入频率f2的值(300<=f2<400f1+f2>600):');
while((f2<300) || (f2>=400))
disp('输入的频率不满足"300<=f2<400"的条件,请重新输入:');
f2=input('f2:');
end
end
disp('请选择求解FFT的算法(输入DIT或者DIF),注意输入字母大写');
select=input('你的选择:''s');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%确定长度N的值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fT=1000;
F1=(300-f1);
F2=(f2-300);
if(F1<F2)
F=F1;
else
F=F2;
end
N1=(fT/F);
i=3;
while ((2^i)<N1)
i=i+1;
end
N=2^(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%得到离散时序列x[n]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(1N);
for n=1:N
x(n)=1.5*cos(2*pi*f1*n/N)+cos(2*pi*300*n/N)+0.5*cos(2*pi*f2*n/N);
end
%画出x[n]的图像
k=1:1:N;
y=1.5*cos(2*pi*f1*k/N)+cos(2*pi*300*k/N)+0.5*cos(2*pi*f2*k/N);
subplot(211);stem(ky'filled');xlabel('n');ylabel('x[n]');title('x[n]的采样图');
m=log2(N);
if(strcmp(select'DIT'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%基2的DIT-FFT算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
change=bin2dec(fliplr(dec2bin([1:N]-1m)))+1;
Y1=x(change);
for s=1:m
Nr=2^s;u=1;
WN=exp(-1j*2*pi/Nr);
for j=1:Nr/2
for k=j:Nr:N
kp=k+Nr/2;
g=Y1(kp)*u;
Y1(kp)=Y1(k)-g;
Y1(k)=Y1(k)+g;
end
u=u*WN;
end
end
Y=Y1;
elseif(strcmp(select'DFT'))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%基2的DIF-FFT算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y2=x;
WN=exp(-1j*2*pi/N);
for s=1:m
for j=1:2^(m-s+1):N-1
for k=1:2^(m-s)
p=j+k;
q=p+2^(m-s);
r=2^(s-1)*mod(p2^(m-s+1));
%g=Y(kp)*u;)
Y2(p)=Y2(p)+Y2(q);
Y2(q)=(Y2(p)+Y2(q))*WN^r;
end
u=u*WN;
end
end
NI=N/2;
for I=1:N-1
if I<NI
t=Y2(I+1);
Y2(I+1)=Y2(NI+1);
Y2(NI+1)=t;
end
T=N/2;
while(NI>=T)
NI=NI-T;
T=T/2;
end
NI=NI+T;
end
% change=bin2dec(fliplr(dec2bin([1:N]-1m)))+1;
% Y=Y2(change);
Y=Y2;
end
subplot(212);stem(changereal(Y)'filled');xlabel('K');ylabel('x[K]');title('x[K]');
disp('是否要继续进行计算?【输入:Y(是)或者N(否)】')
YN=input('是?否?''s');
end
然后进行验证DIF、DIT两种方法的结果是否相同:
- 选择DIT法:
运行结果:
2.选择DIF法:
运行结果:
由上面的验证可知程序的可行性与正确性。
欢迎交流~