快速傅立叶变换(FFT)的计算机实现

信号与系统课程设计——FFT的计算机实现快速傅里叶变换（FFT）的计算机实现赖智鹏华中科技大学电气与电子工程学院0809班U200811806Email:[e-mail]摘要：本文是信号与系统课程的课程设计，旨在熟悉FFT的计算过程，结合DFT物理意义和实验结果加深对傅立叶变换的理解。文章首先用MATLAB对一个简单信号进行FFT仿真，得出频谱图；其次完成了FFT的C语言实现，结合MATLAB作图及数据处理功能得出了C实现下的FFT结果；最后，讨论分析实验结果。关键词：DFT、基--2按时间抽取FFT算法、MATLAB、C、频谱、物理意义1.算法描述1)DFT的运算量1()()()(0,1,2...1)022复数乘：次，复数加：（1）次NnkXkDFTxnxnWkNNnNNNN2)减少运算的方法：化长序列为短序列。如将长度为N的序列分解为两个长度为N/2的序列12()()()(0,1,,1)0222复数乘：次，复数加：（1）次4224NnkXkDFTxnxnWkNNnNNNN利用nkNW的性质（注：本文中的C程序未用到此性质）+rN周期性：=*-共轭对称性：=r可约性：=rnknkWWNNnnWWNNnnWWNN3)C程序采用基--2按时间抽取的FFT算法设输入序列长度为2MN(M为正整数)，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法,也称为Coolkey-Tukey算法。若N不满足条件，则人为地加上若干零值，使'2MN。2.实验设计与步骤为简单起见,同时不失一般性,本实验采用三个余弦成分和一个支流偏置成分叠加所得的信号作为信号源。I.编写Matlab仿真程序2,,,,,,,,,,，其中为被采样信号的01112223330直流成分,,,为第一个正弦成分的参数,,,为第二个正弦成分的参数,,,为第111222333三个正弦成分的参数,为采样频率sII.通过改变2,,,,,,,,,,的输入参数，改变信号波形和采样点0111222333数，观察序列长度FFTaafafaffasafafaffFFTaafafafN或谐波成分对频谱的影响，观察频谱图的变化趋势，并结合频谱图的物理意义作出作出解释解释III.根据Decimation-in-Time(DIT)(Coolkey-Tukey)算法编写C程序，计算结果放在keshe.txt文件输出中，数据调入Matlab中，分析数据,做出频谱图IV.通过改变C程序运行时输入的参数，改变信号波形和采样点数，对比观察频谱的变化趋势，分析输出结果的差异，结合频谱图的物理意义进行解释V.C语言实现下的FFT结果对比MATLAB的仿真结果，观察二者的之间的差异，对之作出合理解释3.MATLAB的FFT仿真和C的FFT实现1)关于频谱的理论分析：Hz,NkN,00N,2特别的，实序列的幅度谱除k=0点外呈轴对称分布Nxfxsxfxskk在1秒内进行采样，频谱分辨率为1注意到频谱的物理意义(在下面的内容中，将利用并一定程度上验证这些性质)I、幅度谱其中为采样点数目直流成分，幅度谱上对应=0的幅度大小为亦即交流成分，幅度谱上对应频率下幅度大小为即II、相位谱直流成分，无相位概念交流成分，相位谱上对应频率下相位的大小就是该频率成分初相位k，。为偶数时，相位谱除k=0N/2点外呈中心对称分布N为奇数时，相位谱除k=0，呈中心对称分布和，；。2)MATLAB的FFT仿真MATLAB程序代码function[output_args]=FFT2(x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs)%这是一个自定义函数，输入被采样的信号的参数和采样频率，然后输出原信号波形、采样信号序列、采样序列幅度谱和相位谱。函数规定信号由一个直流成分（大小为x0），三个余弦成分（各自频率分别为f1,f2,f3,幅值为x1,x2,x3,初相位为w1,w2,w3）,采样频率为fs，采样时间为1秒。x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs作为参数输入t=0:1/fs:1-1/fs;%定义采样时刻N=length(t);%采样序列长度s=x0+x1*cos(2*pi*f1*t+pi/180*w1)+x2*cos(2*pi*f2*t+pi/180*w2)+x3*cos(2*pi*f3*t+pi/180*w3);%采样信号y=fft(s);%快速傅立叶变换tt=0:1/10000:1;%原信号描点ss=x0+x1*cos(2*pi*f1*tt+pi/180*w1)+x2*cos(2*pi*f2*tt+pi/180*w2)+x3*cos(2*pi*f3*tt+pi/180*w3);%原信号figure;plot(tt,ss);grid;title('原始信号');%原信号波形figure;subplot(3,1,1);plot(0:N-1,s(1:N),'-o');xlim([0N-1]);grid;title('采样信号');subplot(3,1,2);plot(0:N-1,abs(y(1:N)),'-o');xlim([0N-1]);grid;xlabel('k');ylabel('幅度');title('理想采样信号的幅度谱');subplot(3,1,3);plot(0:N-1,angle(y(1:N)),'-o');grid;xlabel('k');ylabel('相位');axis([0N-1-pipi]);title('理想采样信号的相位谱');end;改变采样频率依次键入：%其中FFT2()是自定义函数,对信号在1秒内以频率fs进行采样FFT2(1,1,1,90,2,2,180,3,3,180,32)FFT2(1,1,1,90,2,2,180,3,3,180,16)FFT2(1,1,1,90,2,2,180,3,3,180,8)FFT2(1,1,1,90,2,2,180,3,3,180,4)原始信号如图：00.10.20.30.40.50.60.70.80.91-6-4-202468原始信号图1即对信号y=1+cos2/2+2cos43cos6在1秒内分别用32、16、8、4进行采样，然后进行FFT，频谱如下图stttfHzHzHzHz051015202530-10010采样信号051015202530050k幅度理想采样信号序列的幅度谱051015202530-202k幅角理想采样信号序列的相位谱图2、fs=32Hz051015-10010采样信号05101502040k幅度理想采样信号序列的幅度谱051015-202k幅角理想采样信号序列的相位谱图3、fs=16Hz01234567-505采样信号0123456701020k幅度理想采样信号序列的幅度谱01234567-202k幅角理想采样信号序列的相位谱图4、fs=8Hz00.511.522.53-505采样信号00.511.522.530510k幅度理想采样信号序列的幅度谱00.511.522.53-202k幅角理想采样信号序列的相位谱图5、fs=4Hz注意到：1、幅度谱，除却k=0外，图形呈轴对称分布；2、相位谱，N为偶数时除却k=0和N/2外（N为奇数时，除却k=0外），图形呈中心对称分布。理论上,由于复指数的周期性，长度为N（假设N为偶数，N为奇数时类似）的的DFT频谱分析，可把k=N/2,N/2+1,…,N-1看作是负频率-N/2,-N/2+1,-N/2+2,…,-2,-1,特别的，对于实序列而言，正负频率成分的幅值必须相等，而初相位必须相反。下面取fs=8Hz和4Hz部分计算结果分析，其他情况可类似处理。1、8FFT结果的前四个为8.00000000000000+0.00000000000000i，-2.6645352591003e-15+3.99999999999i，-7.9999999999999+1.33226762955e-15i，-12.000000000000+1.065814103640e-14i幅度为8.00000000000000，3.9999999999fHzs9999，7.99999999999999，12.0000000000000此时，N=8，直流成分8.00000000000000/N1,0交流成分3.99999999999999/N/2112,233与=1,=1,=2,=3吻合0123相位为1.57079632679490，3.14159265358979，3.14159265358979化为角度制90.000000xxxxxxxx0000001，180，180.000000000000,即三个频率下初相位分别为90，180，180，与=90,=180,=180吻合1232、4FFT结果的前四个为3.99999999999999+0.00000000000000i，-6.00000000000000+1.99999999999999i，-7.99999999999999+0.00000000000000i，-6.00000000000000-1.99999999999999i幅度为3.99999999999999，6.3245553203fHzs3676，7.99999999999999，6.32455532033676此时，N=4，直流成分3.99999999999999/N=1,0交流成分=6.32455532033676/(N/2)=1.581138830084189,1=1.999999999999998，2=1.581138830084189，3与=1,=1,=2,=3不吻合0123相位为2.819842099xxxxxxxx19315，3.14159265358979，-2.81984209919315化为角度制161.565051177078，180，-161.565051177078，即三个频率下的初相位为161.565051177078，180，-161.565051177078，与=90,=180,=180不吻合123可看到：fs=8Hz时，由频谱结合物理意义计算的得到的余弦成分的幅值和初相位与原信号一致；fs=4Hz时，由频谱结合物理意义计算的得到的余弦成分的幅值和初相位与原信号不同。对此，可根据抽样定理解释。因为f1=1Hz，f2=2Hz，f3=3Hz，所以连续时间信号的截至频率fm=f3=3Hz，所以抽样频率fs=4Hz时，fs2fm=6Hz,所以DTFT频谱出现混叠现象，而DFT实质是对DTFT频谱的频率抽样，所以离散序列的频谱出现混叠现象。直观的，可把混叠现象看作是幅度谱和相位谱中对称的右半部分图形，随着采样频率的减小而往左移动，当采样频率足够小时两部分图形出现重叠交错（结合图6、7、8）。若被抽样信号，截止频率无限大，则混叠现象始终存在。下图从左至右依次为抽样频率fs=8Hz、7Hz、6Hz时的频谱，可看出fs=7Hz为发生混叠的临界频率，即fs=7Hz时恰好不发生混叠现象。01234567-505采样信号0123456701020k幅度理想采样信号序列的幅度谱01234567-202k幅角理想采样信号序列的相位谱图6fs=8Hz不发生混叠0123456-505采样信号012345601020k幅度理想采样信号序列的幅度谱0123456-202k幅角理想采样信号序列的相位谱图7fs=7Hz恰不发生混叠00.511.522.533.544.55-10010采样信号00.511.522.533.544.5501020k幅度理想采样信号序列的幅度谱00.511.522.533.544.55-202k幅角理想采样信号序列的相位谱图8fs=6Hz恰发生混叠改变信号波形键入：%改变信号的直流成分（由1变为2），原