快速傅里叶变换FFT的matlab实现和FFT的简单应用

时，其完成时间为108 （ 为计算机的时钟周期）），故其实现难度是相当大的，同时也严
重制约了 DFT 在信号分析中的应用，故需要提出一种快速的且有效的算法来实现。
正是鉴于 DFT 极其复杂的时间复杂度，1965 年 J.W.Cooley 和 J.W.Tukey 巧妙地利用
WN 因子的周期性和对称性，提出了一个 DFT 的快速算法，即快速傅里叶变换（FFT），从
8
title('时域信号'); 仿真图：
信号与系统课程设计
附件4 源代码：
fs=400;T=1/fs; %采样频率和采样间隔 Tp=0.04;N=Tp*fs; %采样点数 N1=[N,4*N,8*N]; %设定三种截取长度 for m=1:3
n=1:N1(m); xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T); Xk=myfft(xn,4096); fk=[0:4095]/4096/T; subplot(3,2,2*m-1);plot(fk,abs(Xk)/max(abs(Xk))); if m==1 title('矩形窗截取'); end end
N 的 Hamming 窗函数列向量 wn；
（3） 对 x(n) 作 2048 点 FFT 作为 x(t)的近似连续频谱 X ( jf ) .其中 N 为采样点数， N fT ，T 为截取时间长度，取三种长度：0.04s，0,16s,0.32s。
4
信号与系统课程设计
程序源代码及仿真图见附件 4。 四、总结 4.1 对 myfft 实现快速傅里叶变换的评价
5
信号与系统课程设计
六、附件清单
附件 1：
function y=myfft(xr,n) xr=[j*pi*1/8 j*pi*2/8 j*pi*3/8 j*pi*4/8 j*pi*5/8 j*pi*6/8 j*pi*7/8 j*pi*8/8]; n=8; p=0:1:n-1; nu=log2(n); p1=p; b=zeros(1,n); for t=1:nu;
信号与系统课程设计
快速傅里叶变换 FFT 的 matlab 实现和 FFT 的简单应用
（1卿立艳 ２００８１０ቤተ መጻሕፍቲ ባይዱ０１８
０８级 电子二班）
【摘要】 在信号处理中，DFT（离散傅里叶变换）的计算具有举足轻重的地位。但是基于 其复杂的计算，直接应用起来十分麻烦，基于此，本文利用 Matlab 软件对有限长度信号的 DFT 进行改进，提出 FFT（快速傅里叶变换），并利用 FFT 对所给连续时间和离散时间信号 做了频谱分析。
本文基于时间抽选奇偶分解算法对 DFT 进行了改进，提出了快速傅里叶变换 FFT，并 用 Matlab 实现了 FFT，并用其对所给信号进行了频谱分析。在第一个实例里，利用 myfft
函数对离散信号 x[n] 进行了傅里叶变换得到它的频域函数和时域函数，从仿真图可以看出
myfft 实现的傅里叶变换时成功的。在第二个实例里面，也是利用 myfft 对信号序列进行傅 里叶变换，并对信号利用窗函数减少频谱间的干扰，较准确地分析了所给信号的频谱。 4.2 其它 FFT 算法简介
h=2^(m-1); k=1; while(k<n+1) for t=1:h;
y=bitshift(k-1,nu-m,nu)+1; %ÇówµÄÃÝ´ÎÊý xch(k)=xr(k)+w(y)*xr(k+h); k=k+1; end; for t=1:h;
y=bitshift(k-1-h,nu-m,nu)+1; %ÇówµÄÃÝ´ÎÊý xch(k)=xr(k-h)-xr(k)*w(y); k=k+1; end; end; xr=xch; % ¼ÆËãx(k)½áÊø end; y=xr
p2=floor(p1/2); b=b*2+(p1-2*p2); p1=p2; end; yr(p+1)=xr(b+1); xr=yr; % 倒位序结束 t=0:n/2-1; %计算因子 w 开始 (只计算 w0 到 w n/2-1) for v=0:n/2-1; w=exp(-2*i*pi*t/n); end; %计算因子 w 结束 for m=1:nu;% 计算 x(k)开始 h=2^(m-1); k=1; while(k<n+1)
设 x(t) cos(200t) sin(100t) cos(50t) ，利用傅里叶变换分析其频谱结构，选 择不同的截取长度 T ，观察存在的截取效应，并试用加窗的方法减少谱间干扰。
选取的参数为： （1） 频率 f 400Hz,T 1/ f ；
（2） 采样信号序列 x[n] x(nT )w(n) ， w(n) 是窗函数，选取两种窗函数：矩形窗函数 w(n) RN (n) 和 Hamming 窗，后者在程序中调用函数 w(n)=hamming(N)产生成都为
得到结果：
6
信号与系统课程设计
附件 2 一个验证 myfft 正确与否的小程序： 代码： N=8; n=0:N-1; xn=cos(pi*n/8); Xk=myfft(xn,N); plot(Xk);stem(n,abs(Xk),'.');axis([0,20,0,20]);ylabel('|Xk|'); title('8点FFT变换'); 仿真图：
7
信号与系统课程设计
附件 3:源代码： %产生两个正弦加白噪声 N=2^8; f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N); %应用FFT求频谱 subplot(2,2,1); plot(x(1:N/4)); title('原始信号'); f=-0.5:1/N:0.5-1/N; x=myfft(x); y=ifft(x); subplot(2,2,2); plot(f,fftshift(abs(x))); title('频域信号'); subplot(2,2,3); plot(real(x(1:N/4)));
本文其实还只是快速傅里叶变换的一个简单实现而已。快速傅里叶变换，顾名思义， 它应该还有其它实现的算法。本文用的是基于时间抽取奇偶分解的算法，在常规的快速傅里 叶变换实现中，还有其它两种非常重要的算法：基于频率抽选奇偶分解算法和混合基算法。 其中混合基算法是运用最广泛的一种 FFT 算法，并在很多领域取得了广泛的应用。由于笔者 水平有限，便不再对上述两种算法进行深入的讨论，留待以后慢慢尝试！ 4.3 FFT 的进一步深入应用
N 1
X (k) x[n]WNnk n0
其中 WNnk
j 2 nk
e N 。
对任意 0 m N 1，
N 1
X (m) x[n]WNnm x[0]WN0m x[1]WN1m ... x[N 1]WN(N1)m n0
1
（1） （2）
信号与系统课程设计
nu=log2(n); p1=p; b=zeros(1,n); for t=1:nu;
只能用一个字来形容：无奈，累！本来期末时间都要复习考试，可是还要做这个课程设 计，花掉我们很多的时间和精力，最后还搞出来一个不是十分满意的东西。哎，一个字：苦 不堪言！
不过，也并非只是感觉苦不堪言，也有感觉高兴的地方：学了 Matlab 了，小有收获， 还熬了夜，提前体验了一把工作加班到深夜的滋味，很不爽，但又很爽！
五、参考文献： [1] 余成波，陶红艳。数字信号处理及 MATLAB 实现，北京：清华大学出版社，2008 [2]（美）Edward W.Kamen, Bonnie S.Heck 著，高强译。 信号与系统基础教程，北京：电
子工业出版社，2007 [3] 曹弋，赵阳。MATLAB 实用教程，北京：电子工业出版社，2007
3
信号与系统课程设计
end; xr=xch; end;% 计算 x(k)结束 y=xr%输出变换后的结果 程序结束 附件清单 1 和 2 验证了其正确性。
三、利用快速傅里叶变换 myfft 函数实现频域分析实例：
3.1 设 x[n]是由两个正弦信号即白噪声的叠加，请用傅里叶变换对其作频域分析。
x[n] 的产生：
关 键 词：DFT，FFT，有限长度信号，频谱分析。
一、前言： 傅里叶变换在信号处理中具有十分重要的作用，但是基于离散时间的傅里叶变换具有很
大的时间复杂度，根据傅里叶变换理论，对一个有限长度且长度为 N 的离散信号，做傅里 叶变换的时间复杂度为 O(N 2 ) ，当 N 很大时，其实现的时间是相当惊人的（比如当 N 为104
N=2^8; f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N); 其傅里叶变换和频域分析，以及仿真图详见附件 3。
3.2 利用快速傅里叶变换 FFT 对连续信号作谱分析。
由于傅里叶变换在信号处理中具有举足轻重的作用，故 FFT 有着非常广泛的应用。下 一步的工作便是利用 myfft 算法设计一个 fir 滤波器，多所给参杂了噪声的语音信号进行滤 波，由于这个设计需要用到其它数字信号处理的知识，可能要等到下学期学习了数字信号课 以后才能实现了。 4.4 课程设计之后的感悟
仿真图：
10
for t=1:h; y=bitshift(k-1,nu-m,nu)+1; %求 w 的幂次数 xch(k)=xr(k)+w(y)*xr(k+h); k=k+1;
end; for t=1:h;