对于FFT和IFFT算法和频谱分析研究报告报告

1、-. z.对于FFT和IFFT的算法和频谱分析的研究（The algorithms and spectrum analysis ofFFT and IFFT）摘要：目的在于研究前人的工作结果，对FFT和IFFT有更清楚的认识。主要通过MATLAB的编程完成对FFT和IFFT的算法和频谱分析。首先通过matlab的编程实现FFT和IFFT的这两个函数。然后用已经编译成功的函数实现升余弦滚降。用FFT分析三角函数和三角波函数。用IFFT将上述结果重新变回到时域，通过作图分析变换前后信号的差异。得出了关于fft和ifft函数的分析和关于三角函数和三角波函数的频谱分析的结论关键词：MATLAB FFT

2、 IFFT 升余弦滚降函数 三角函数 三角波函数Abstract：pleted the main algorithm and spectral analysis of FFT and IFFT by MATLAB programming. First, through the MATLAB programming to achieve the two functions FFT and IFFT. Then use has been successfully piled function raised cosine. Analysis of trigonometric function and

3、 triangle function by FFT. With IFFT the results back in time domain, by mapping differences before and after signal transformation.Key words: MATLAB ,FFT ,IFFT,Raised cosine function,Trigonometric,Triangular wave function引言1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在计算数学杂志上发表了机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的

4、一种快速有效的计算方法的文章。它的思路建立在对DFT运算在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（FastFourierTransform），简称FFT，1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法， 2使运算效率进一步提高。FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。随

5、着科学的进步，FFT算法的重要意义已经远远超过傅里叶分析本身的应用。FFT算法之所以快速，其根本原因在于原始变化矩阵的多余行，此特性也适用于傅里叶变换外的其他一些正交变换，例如，快速沃尔什变换、数论变换等等。在FFT的影响下，人们对于广义的快速正交变换进行了深入研究，使各种快速变换在数字信号处理中占据了重要地位。因此说FFT对数字信号处理技术的发展起了重大推动作用。快速傅里叶变换(FastFourier Tranformation，FFT)是将一个大点数N的DFT分解为若干小点的DFT的组合。将运算工作量明显降低，从而大大提高离散傅里叶变换(DFT)的计算速度,从而更加适合进行实时运算。因各个

6、科学技术领域广泛的使用了FFT技术它大大推动了信号处理技术的进步，现已成为数字信号处理强有力的工具，本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点，并完成了基于MATLAB的实现。最后通过FFT和IFFT的两个应用升余弦滚降和确定函数的频谱分析来分别验证FFT和IFFT的正确性和优越性。FFT的算法 1 1.1 FFT算法的基本思想设离散的有限长时间序列*(n)，0nN-1，则其离散傅立叶变换为：这样，矩阵W 中有许多相同的元素，从而可以简化DFT的运算过程FFT 算法有许多形式，笔者只讨论最基本的时间抽取基-2FFT 算法1.2 算法分析一个N 点长序列，直接用DFT 方法需要复数乘

7、法N次；复数加法N(N-1)次。而由图2 可知，采用FFT 则只需要复数乘法次；复数加法 次。当 时， 这样，运算速度提高了1-2 个数量级图1 为FFT 算法和直接DFT 算法所需运算量与计算点数N 的关系曲线显然，N 越大时，优越性越明显但当N 相当大时，利用单机串行进行FFT 运算同样满足不了实时系统的需要 1 1.3 算法的程序实现思想及分析首先检验待变换的序列的元素个数是否为2的幂次方个，如果不是的话则将其补零使之成为2的整数幂次方个。然后利用已经编好的位倒序子程序输出位倒序序号，将输入序列不断分组，进行处理后从新分组，直至完成最后的处理即可输出变换后的结果。需要注意的是，这里fft

8、的变换后结果的元素个数可能与原输入序列的个数不一样，因为如果不是2的幂次方个的话输入序列后面是要补零的。（主程序为fft_dit和fft_dif）1.4流程示意图整个FFT 频谱分析与显示过程可用图2 所示的流程图示意.fft按频域抽取算法fft按时域抽取算法图22.IFFT的频谱变换基本原理在实验中为了简化算法我们直接利用前面已经编好的fft_dit或fft_dif这两个现成函数来实现，基本原理如下：将要变换的频谱序列先取共轭然后将其送入前面的函数，将变化后的结果再取共轭即实现IFFT的功能。主程序为ifft_my图33升余弦滚降3.1升余弦滚降原理要实现无码间干扰基带传输时，系统必须满足奈

9、奎斯特准则即：对于上述公式，我们分3种情况来说明其含义：Ts1/2W情况，Z（f）由频率间隔为1/Ts的*（f）曲线无频率重叠地周期性复制并相加构成的，它还是周期性频谱。在这种情况下，有一特定频谱可满足无码间干扰传输的条件，它就是已获广泛应用的升余弦谱。升余弦滤波器的传递函数表示式为：称为滚降因子，取值为01。在=0时，滤波器的带宽W为1/(2Ts),称为奈奎斯特带宽；=0.5时，滤波器的截止频率W=（1+）/（2Ts）=0.75Rs;=1时，滤波器的截止频率W=Rs。 3 43.2 升余弦滚降的实现及其与快速傅里叶逆变换的关系升余弦滚降函数是在基带无码间干扰传输中经常用到的频域函数。其主要特

10、性是升余弦滚降函数经过频域平移叠加后能够成一个在各个频域幅度恒定的频域函数。我们在实现升余弦滚降时可以首先在频域实现=0，=0.5，=0.75和=1四种升余弦滚降函数，然后通过自己编写的ifft_my进行傅里叶逆变换，并将变换后的时域结果反映在图上，分别对比不同是对应的时域上的不同的时域特性。3.3 变换前的时域特性和变换后的频域特性 图44 利用fft和ifft进行具体函数频谱分析的实例4.1 三角函数的频谱分析及其信号的恢复在实际信号的分析中，三角函数是非常常见和基本的信号，在这里我们就对三角函数进行分析。在分析中我们会碰到要分析的函数的采样点数不是2的整数次幂个，我们要对它进行补零处理。

11、4.2 三角波函数的分析在分析中因为可能会用到自己编写的T2F函数和F2T函数，但是这两个函数仍然是建立在自己编写的fft_dit和ifft_my的基础上的。只是对输入变量进行了更加完善的处理，一个是补零，另一个是对频域进行了搬移，将原来pi2pi的部分搬到-pi0，因为大家在看频谱时比较习惯频谱是关于零对称的。4.3 分析结果5 结论5.1 关于fft和ifft函数的分析在fft和ifft的实现过程中，的确能降低运算的次数，但是也正好印证了一个很著名的理论，时间换空间，空间换时间”。在fft和ifft的算法中，我们降低运算的次数是以占用更多的空间换来的。每次将要变换的的序列进行分组，然后对每

12、个组进行处理，虽然降低了运算次数，但是也增加了运算空间的占用。5.2 关于升余弦滚降函数的结论通过图形我们可以观察出，对于不同的函数，时域主要部分占用的宽度是一样的。但是随着的增大，时域的起伏是越来越小的。因此我们可以认为，越大，时域起伏越小，因此在非线性系统中传输时的失真越小。但是与此同时带来的缺点是占用的越大，频域占用的带宽是越大的，越利于在限带系统中传输，频率的利用率也就越低。5.3 关于三角函数和三角波函数的频谱分析结论（1）根据理论分析，三角函数的频谱因该是几个冲激函数的组合，三角波函数的频谱应该是Sa函数平方的形式。而在实际分析中我们通过观察频谱图可以看出分析的结果基本符合理论分析

13、结果。这也侧面印证了我们自己编写的fft_dit和ifft_my的正确性和有效性。（2）在通过频谱恢复时域时，我们观察到恢复后的三角函数的时域函数有所失真，三角波函数基本无失真。原因可能是在T2F和F2T函数中可能因为补零改变了元素个数，因此引起了恢复后的函数时域出现失真。致：谷群 若寰 王敏附录1： 1 冬初,何 飞.FFT算法的并行处理研究.城市学院学报(自然科学版) .2005-06-30 2Oppenheim A V, Schafer R W.Digit al Signal Proces sing M . N Y: Prent ice-Hall , 1975.3Proak is J G

14、, Manolakis D G. Digital Signal Proces sing-PrinciplesM . Algorithms and applicat ions ( 2nd Edit ion) , N Y:Print ice-Hall, 1995. 4灵智. M atlab 在数字信号处理课程设计中的应用 J. 水利职业学院刊, 2008. 3: 11- 12.附录2：主要程序和代码升余弦滚降：Fftfunction num3=fft_dif(in) N=length(in); num=in; num3=zeros(1,N); k=log2(N); for n=1:k for c=

15、1:2(n-1) num1=zeros(1,2(k+1-n); num2=zeros(1,2(k+1-n); num1=num(c-1)*2(k+1-n)+1:c*2(k+1-n); for m=1:2(k-n) num2(m)=num1(m)+num1(m+2(k-n) num2(m+2(k-n)=(num1(m)-num1(m+2(k-n)*e*p(-i*(m-1)*2*pi/2(k+1-n); end num(c-1)*2(k+1-n)+1:c*2(k+1-n)=num2; end end for n=1:N num3(n)=num(weidao*u(n,N)+1); endfuncti

16、on num3=fft_dit(in) N=length(in); num=zeros(1,N); num3=zeros(1,N); k=log2(N); for n=1:N num(n)=in(weidao*u(n,N)+1); end for n=1:k for c=1:2(k-n) num1=zeros(1,2n); num2=zeros(1,2n); num1=num(c-1)*2n+1:c*2n); for m=1:2(n-1) num2(m)=num1(m)+num1(m+2(n-1)*e*p(-i*(m-1)*2*pi/2n); num2(m+2(n-1)=num1(m)-num

17、1(m+2(n-1)*e*p(-i*(m-1)*2*pi/2n); end num(c-1)*2n+1:c*2n)=num2; end end %num3(1)=num(N); %num3(2:N)=num(1:N-1); num3=num;IFFTfunction out=ifft_my(in) num1=conj(in); num2=fft_dif(num1);out=conj(num2)/length(in);%通过FFT和IFFT来分析升余弦滚降函数的时域和频域的特点clear all;a=0;N=128;l=2;figure(1)f,out_f=shengyu*iangunjiang

18、(a,N,l);subplot(221);plot(f,out_f);*label(f);ylabel(a=0时的频谱幅度);a*is(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(222);plot(f,out_t_);*label(t);ylabel(a=0时的时域幅度);a*is(-0.25 0.25 -0.1 0.5);a=0.5;f,out_f=shengyu*iangunjiang(a,N,l);sub

19、plot(223);plot(f,out_f);*label(f);ylabel(a=0.5时的频谱幅度);a*is(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(224);plot(f,out_t_);*label(t);ylabel(a=0.5时的时域幅度);a*is(-0.25 0.25 -0.1 0.5);a=0.75;figure(2)f,out_f=shengyu*iangunjiang(a,N,l);

20、subplot(221);plot(f,out_f);*label(f);ylabel(a=0.75时的频谱幅度);a*is(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(222);plot(f,out_t_);*label(t);ylabel(a=0.75时的时域幅度);a*is(-0.25 0.25 -0.1 0.5);a=1;f,out_f=shengyu*iangunjiang(a,N,l);subplot

21、(223);plot(f,out_f);*label(f);ylabel(a=1时的频谱幅度);a*is(-2 2 0 1.5);out_f_=out_f(N/2+1):N),out_f(1:N/2);out_t=ifft_my(out_f_);out_t_=out_t(N/2+1):N),out_t(1:N/2);subplot(224);plot(f,out_t_);*label(t);ylabel(a=1时的时域幅度);a*is(-0.25 0.25 -0.1 0.5);逆位倒序：function out=niweidao*u(in,long) l=log2(long); num1=ze

22、ros(1,l); %out=0; for i=1:l num1(i)=mod(in,2); out=(in-mod(in,2)/2; end % for i=1:l % num2(i)=num1(l+1-i); % end % for i=1:l % out=(2(i-1)*num1(i)+out; % endRAND01： function s=rand01(p,m,n) % 输入参数: % p:0-1分布中1的概率 % m,n:产生的随机变量样本个数mn% 输出:产生的随机变量样本矢量升余弦滚降：*=rand(m,n);s=(sign(*-p+eps)+1)/2;% eps = 2(-5

23、2).function f,out_f=shengyu*iangunjiang(a,N,l) f=-l:2*l/N:l-2*l/N; out_f=zeros(1,N); out_f(1:N/(2*l)*(1-a)=0; k=-pi/2:l*pi/(a*N):pi/2-l*pi/(a*N) out_f(N/(2*l)*(1-a)+1:N/(2*l)*(1+a)=(sin(k)+1)/2; out_f(N/(2*l)*(1+a)+1:N/2)=1; for m=1:N/2 out_f(N+1-m)=out_f(m); end位倒叙function out=weidao*u(in1,long) in

24、=in1-1; l=log2(long); num1=zeros(1,l); out=0; for i=1:l num1(i)=mod(in,2); in=(in-mod(in,2)/2; end for i=1:l num2(i)=num1(l+1-i); end for i=1:l out=(2(i-1)*num2(i)+out;end自余弦滚降function Ra=zi*iangguan(in) N=length(in); Ra=zeros(1,N); for m=1:N for k=1:N-1-m Ra(m)=in(k)*in(k+m-1)+Ra(m); end Ra(m)=Ra(m

25、)/(N+1-m); end频谱分析F变时域function m=F2T(M,fs)%-输入参数%M：信号的频谱%fs:系统采样频率%-输出(返回)参数%m:傅里叶逆变换后的信号，注意其长度为2的整数次幂，利用其画波形时，要注意选取m的一部分，选取长度和所给时间序列t的长度要一致，plot(t,m(1:length(t)，否则会出错。m = real(ifft(M)*fs;FFTfunction num3=fft_dif(in1,k)if nargin=1 N=length(in1);else N=k;end in=zeros(1,N); in(1:length(in1)=in1; num=i

26、n; num3=zeros(1,N); k=log2(N); for n=1:k for c=1:2(n-1) num1=zeros(1,2(k+1-n); num2=zeros(1,2(k+1-n); num1=num(c-1)*2(k+1-n)+1:c*2(k+1-n); for m=1:2(k-n) num2(m)=num1(m)+num1(m+2(k-n) num2(m+2(k-n)=(num1(m)-num1(m+2(k-n)*e*p(-i*(m-1)*2*pi/2(k+1-n); end num(c-1)*2(k+1-n)+1:c*2(k+1-n)=num2; end end fo

27、r n=1:N num3(n)=num(weidao*u(n,N)+1); endfunction num3=fft_dit(in1,k)if nargin=1 N=length(in1);else N=k;end in=zeros(1,N); in(1:length(in1)=in1; num=zeros(1,N); num3=zeros(1,N); k=log2(N); for n=1:N num(n)=in(weidao*u(n,N)+1); end for n=1:k for c=1:2(k-n) num1=zeros(1,2n); num2=zeros(1,2n); num1=num

28、(c-1)*2n+1:c*2n); for m=1:2(n-1) num2(m)=num1(m)+num1(m+2(n-1)*e*p(-i*(m-1)*2*pi/2n); num2(m+2(n-1)=num1(m)-num1(m+2(n-1)*e*p(-i*(m-1)*2*pi/2n); end num(c-1)*2n+1:c*2n)=num2; end end %num3(1)=num(N); %num3(2:N)=num(1:N-1); num3=num;function M,m,df=fftseq(m,ts,df)%各参数含义与子函数T2F中的完全相同，完成fs = 1/ts;if na

29、rgin =2 n1 =0;else n1 = fs/df;end n2 = length(m);n = 2(ma*(ne*tpow2(n1),ne*tpow2(n2);M = fft_dit(m,n);M=M(n/2+1:n),M(1:n/2);m = m,zeros(1,n-n2);df = fs/n;function out=ifft_my(in) num1=conj(in); num2=fft_dit(num1); out=conj(num2)/length(in);function out=niweidao*u(in,long) l=log2(long); num1=zeros(1,

30、l); %out=0; for i=1:l num1(i)=mod(in,2); out=(in-mod(in,2)/2; end % for i=1:l % num2(i)=num1(l+1-i); % end % for i=1:l % out=(2(i-1)*num1(i)+out; % endclear alldt=0.001;df=0.2;fs=1/dt;t=-10:dt:10-dt;y=sin(2*pi*100*t)+2*sin(2*pi*400*t);M1,m1,df1,f1=T2F(y,dt,df,fs);figure(1)subplot(211);plot(t,y);*label(t);ylabel(三角函数频谱分析前的时域);a*is(-0.05,0.05,-5,5);subplot(212);plot(f1,abs(M1);*label(f);ylabel(三角函数频谱分析后的频谱);figure(2)z=zeros(1,length(t);z=zeros(1,length(t)/4),0:5*4/length(t):5-5*4/length(t),zeros(1,length(t)/2);for n=1:length(t)/2 z(length(t)+1-