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

《对于FFT和IFFT的算法和频谱分析的研究》由会员分享，可在线阅读，更多相关《对于FFT和IFFT的算法和频谱分析的研究（25页珍藏版）》请在金锄头文库上搜索。

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

2、的分析和关于三角函数和三角波函数的频谱分析的结论 关键词：MATLAB FFT IFFT 升余弦滚降函数 三角函数 三角波函数 Abstract： Completed 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 compiled function raised

3、 cosine. Analysis of trigonometric function and 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

4、.W.Tukey）在计算数学杂志上发表了“机器计算傅立叶级数的一种算法”的文章， 这是一篇关于计算 DFT的一种快速有效的计算方法的文章。它的思路建立在对 DFT运算内在规律的认识之上。这篇文章的发表使 DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（Fast Fourier Transform） ，简称 FFT，1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法， 2使运算效率进一步提高。FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改

5、进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。随着科学的进步， FFT算法的重要意义已经远远超过傅里叶分析本身的应用。FFT算法之所以快速，其根本原因在于原始变化矩阵的多余行， 此特性也适用于傅里叶变换外的其他一些正交变换，例如，快速沃尔什变换、数论变换等等。在 FFT的影响下，人们对于广义的快速正交变换进行了深入研究， 使各种快速变换在数字信号处理中占据了重要地位。因此说 FFT对数字信号处理技术的发展起了重大推动作用。 快速傅里叶变换(FastFourier Tranformation，FFT)是将一个大点数 N的

6、DFT分解为若干小点的 DFT的组合。 将运算工作量明显降低， 从而大大提高离散傅里叶变换(DFT)的计算速度, 从而更加适合进行实时运算。 因各个科学技术领域广泛的使用了 FFT技术它大大推动了信号处理技术的进步， 现已成为数字信号处理强对于 FFT和 IFFT的算法和频谱分析的研究 有力的工具，本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点，并完成了基于 MATLAB的实现。 最后通过 FFT和 IFFT的两个应用升余弦滚降和确定函数的频谱分析来分别验证 FFT和 IFFT的正确性和优越性。 1.FFT的算法 1 1.1 FFT算法的基本思想 设离散的有限长时间序列x(n)， 0

7、n N-1，则其离散傅立叶变换为： 这样，矩阵W 中有许多相同的元素，从而可 以简化DFT的运算过程FFT 算法有许多形式， 笔者只讨论最基本的时间抽取基-2FFT 算法 1.2 算法分析 一个N 点长序列，直接用DFT 方法需要复数乘法N 次；复数加法N(N-1)次。而由图2 可知，采用FFT 则只需要复数乘法次；复数加法 次。当 时， 这样，运算速度提高了1-2 个数量级图1 为FFT 算法和直接DFT 算 法所需运算量与计算点数N 的关系曲线显然，N 越大时，优越性越明显但当N 相当大时，利用单机串行进行 FFT 运算同样满足不了实时系统的需要 1 对于 FFT和 IFFT的算法和频谱分

8、析的研究 1.3 算法的程序实现思想及分析 首先检验待变换的序列的元素个数是否为 2 的幂次方个，如果不是的话则将其补零使之成为 2 的整数幂次方个。 然后利用已经编好的位倒序子程序输出位倒序序号，将输入序列不断分组，进行处理后从新分组，直至完成最后的处理即可输出变换后的结果。需要注意的是，这里 fft的变换后结果的元素个数可能与原输入序列的个数不一样， 因为如果不是 2 的幂次方个的话输入序列后面是要补零的。 （主程序为 fft_dit和 fft_dif） 1.4流程示意图 整个FFT 频谱分析与显示过程可用图2 所示的流程图示意. x(0)X(0)x(1)x(2)x(3)x(4)x(5)x

9、(6)x(7)X(5)X(4)X(1)X(2)X(6)X(3)X(7) fft按频域抽取算法 对于 FFT和 IFFT的算法和频谱分析的研究 x(0)x(4)x(6)x(2)x(1)x(5)x(7)x(3)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7) fft按时域抽取算法 图2 2.IFFT的频谱变换基本原理 在实验中为了简化算法我们直接利用前面已经编好的fft_dit或fft_dif这两个现成函数来实现， 基本原理如下： 将要变换的频谱序列先取共轭然后将其送入前面的函数，将变化后的结果再取共轭即实现 IFFT的功能。主程序为ifft_my 待变换序列取共轭快速傅里叶变换

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

11、 （2Ts）=0.75Rs; =1时，滤波器的截止频率W=Rs。 3 4 3.2 升余弦滚降的实现及其与快速傅里叶逆变换的关系 升余弦滚降函数是在基带无码间干扰传输中经常用到的频域函数。其主要特性是升余弦滚降函数经过频域平移叠加后能够成一个在各个频域幅度恒定的频域函数。我们在实现升余弦滚降时可以首先在频域实现 =0， =0.5， =0.75和 =1四种升余弦滚降函数，然后通过自己编写的ifft_my进行傅里叶逆变换，并将变换后的时域结果反映在图上， 分别对比 不同是对应的时域上的不同的时域特性。 3.3 变换前的时域特性和变换后的频域特性 对于 FFT和 IFFT的算法和频谱分析的研究 图4

12、4 利用fft和ifft进行具体函数频谱分析的实例 4.1 三角函数的频谱分析及其信号的恢复 在实际信号的分析中，三角函数是非常常见和基本的信号，在这里我们就对三角函数进行分析。在分析中我们会碰到要分析的函数的采样点数不是2 的整数次幂个，我们要对它进行补零处理。 4.2 三角波函数的分析 在分析中因为可能会用到自己编写的T2F函数和F2T函数，但是这两个函数仍然是建立在自己编写的fft_dit和ifft_my的基础上的。 只是对输入变量进行了更加完善的处理，一个是补零，另一个是对频域进行了搬移，将原来pi2pi的部分对于 FFT和 IFFT的算法和频谱分析的研究 搬到-pi0，因为大家在看频

13、谱时比较习惯频谱是关于零对称的。 4.3 分析结果 对于 FFT和 IFFT的算法和频谱分析的研究 5 结论 5.1 关于fft和ifft函数的分析 在fft和ifft的实现过程中，的确能降低运算的次数，但是也正好印证了一个很著名的理论， “时间换空间，空间换时间” 。在fft和ifft的算法中，我们降低运算的次数是以占用更多的空间换来的。 每次将要变换的的序列进行分组， 然后对每个组进行处理，虽然降低了运算次数，但是也增加了运算空间的占用。 5.2 关于升余弦滚降函数的结论 通过图形我们可以观察出，对于不同 的函数，时域主要部分占用的宽度是一样的。但是随着 的增大，时域的起伏是越来越小的。因

14、此我们可以认为，越大，时域起伏越小，因此在非线性系统中传输时的失真越小。但是与此同时带来的缺点是占用的 越大， 频域占用的带宽是越大的， 越利于在限带系统中传输，频率的利用率也就越低。 5.3 关于三角函数和三角波函数的频谱分析结论 （1 ）根据理论分析，三角函数的频谱因该是几个冲激函数的组合，三角波函数的频谱应该是Sa函数平方的形式。 而在实际分析中我们通过观察频谱图可以看出分析的结果基本符合理论分析结果。这也侧面印证了我们自己编写的fft_dit和ifft_my的正确性和有效性。 （2 ）在通过频谱恢复时域时，我们观察到恢复后的三角函数的时域函数有所失真，三角波函数基本无失真。原因可能是在

15、 T2F和F2T 函数中可能因为补零改变了元素个数，因此引起了恢复后的函数时域出现失真。 致谢： 谷群 陈若寰 王敏 对于 FFT和 IFFT的算法和频谱分析的研究 附录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, Manolakis D G. Digital Signal Proces sing-PrinciplesM .

16、 Algorithms and applicat ions ( 2nd Edit ion) , N Y:Print ice-Hall, 1995. 4崔灵智. M atlab 在数字信号处理课程设计中的应用 J. 山东水利职业学院刊, 2008. 3: 11- 12. 附录2：主要程序和代码 升余弦滚降： Fft function num3=fft_dif(in) N=length(in); num=in; 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

17、-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)*exp(-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(weidaoxu(n,N)+1); end function num3=fft_dit(in) N=length(in); num=zeros(1,N);

18、 num3=zeros(1,N); k=log2(N); for n=1:N num(n)=in(weidaoxu(n,N)+1); end for n=1:k for c=1:2(k-n) num1=zeros(1,2n); num2=zeros(1,2n); 对于 FFT和 IFFT的算法和频谱分析的研究 num1=num(c-1)*2n+1:c*2n); for m=1:2(n-1) num2(m)=num1(m)+num1(m+2(n-1)*exp(-i*(m-1)*2*pi/2n); num2(m+2(n-1)=num1(m)-num1(m+2(n-1)*exp(-i*(m-1)*2

19、*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; IFFT function 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=shengyuxiangunjiang(a,N,l); su

20、bplot(221); plot(f,out_f); xlabel(f); ylabel(a=0 时的频谱幅度); 对于 FFT和 IFFT的算法和频谱分析的研究 axis(-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_); xlabel(t); ylabel(a=0 时的时域幅度); axis(-0.25 0.25 -0.1 0.5); a=0.5; f,out_f=

21、shengyuxiangunjiang(a,N,l); subplot(223); plot(f,out_f); xlabel(f); ylabel(a=0.5 时的频谱幅度); axis(-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); 对于 FFT和 IFFT的算法和频谱分析的研究 subplot(224); plot(f,out_t_); xlabel(t); ylabel(a=0.5 时的时域幅度); axis(-0.

22、25 0.25 -0.1 0.5); a=0.75; figure(2) f,out_f=shengyuxiangunjiang(a,N,l); subplot(221); plot(f,out_f); xlabel(f); ylabel(a=0.75 时的频谱幅度); axis(-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_); 对于 FFT和 IFFT的算法和频谱分

23、析的研究 xlabel(t); ylabel(a=0.75 时的时域幅度); axis(-0.25 0.25 -0.1 0.5); a=1; f,out_f=shengyuxiangunjiang(a,N,l); subplot(223); plot(f,out_f); xlabel(f); ylabel(a=1 时的频谱幅度); axis(-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

24、(f,out_t_); xlabel(t); ylabel(a=1 时的时域幅度); axis(-0.25 0.25 -0.1 0.5); 逆位倒序： function out=niweidaoxu(in,long) 对于 FFT和 IFFT的算法和频谱分析的研究 l=log2(long); num1=zeros(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)

25、+out; % end RAND01： function s=rand01(p,m,n) % 输入参数: % p:0-1 分布中 1 的概率 % m,n: 产生的随机变量样本个数 m n % 输出:产生的随机变量样本矢量 升余弦滚降： x=rand(m,n); s=(sign(x-p+eps)+1)/2; % eps = 2(-52). 对于 FFT和 IFFT的算法和频谱分析的研究 function f,out_f=shengyuxiangunjiang(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;

26、 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=weidaoxu(in1,long) in=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:

27、l num2(i)=num1(l+1-i); 对于 FFT和 IFFT的算法和频谱分析的研究 end for i=1:l out=(2(i-1)*num2(i)+out; end 自余弦滚降 function Ra=zixiangguan(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)/(N+1-m); end 频谱分析 F变时域 function m=F2T(M,fs) %-输入参数 %M：信号的频谱 %fs:系统采样频率 %-输出(返回

28、)参数 对于 FFT和 IFFT的算法和频谱分析的研究 %m:傅里叶逆变换后的信号， 注意其长度为 2 的整数次幂， 利用其画波形时， 要注意选取 m 的一部分，选取长度和所给时间序列 t 的长度要一致，plot(t,m(1:length(t)，否则会出错。 m = real(ifft(M)*fs; FFT function 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=in; num3=zeros(1,N); k=log2(N); fo

29、r 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); 对于 FFT和 IFFT的算法和频谱分析的研究 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)*exp(-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 n

30、um3(n)=num(weidaoxu(n,N)+1); end function 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); 对于 FFT和 IFFT的算法和频谱分析的研究 k=log2(N); for n=1:N num(n)=in(weidaoxu(n,N)+1); end for n=1:k for c=1:2(k-n) num1=zeros(1,2n); num2=ze

31、ros(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)*exp(-i*(m-1)*2*pi/2n); num2(m+2(n-1)=num1(m)-num1(m+2(n-1)*exp(-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 中

32、的完全相同，完成 fs = 1/ts; 对于 FFT和 IFFT的算法和频谱分析的研究 if nargin =2 n1 =0; else n1 = fs/df; end n2 = length(m); n = 2(max(nextpow2(n1),nextpow2(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); fu

33、nction out=niweidaoxu(in,long) l=log2(long); num1=zeros(1,l); %out=0; for i=1:l 对于 FFT和 IFFT的算法和频谱分析的研究 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; % end clear all dt=0.001; df=0.2; fs=1/dt; t=-10:dt:10-dt; y=sin(2*pi*10

34、0*t)+2*sin(2*pi*400*t); M1,m1,df1,f1=T2F(y,dt,df,fs); figure(1) subplot(211); plot(t,y); xlabel(t); ylabel(三角函数频谱分析前的时域); 对于 FFT和 IFFT的算法和频谱分析的研究 axis(-0.05,0.05,-5,5); subplot(212); plot(f1,abs(M1); xlabel(f); ylabel(三角函数频谱分析后的频谱); figure(2) z=zeros(1,length(t); z=zeros(1,length(t)/4),0:5*4/length(

35、t):5-5*4/length(t),zeros(1,length(t)/2); for n=1:length(t)/2 z(length(t)+1-n)=z(n); end M2,m2,df2,f2=T2F(z,dt,df,fs); subplot(211); plot(t,z); xlabel(t); ylabel(三角波频谱分析前的时域); axis(-10,10,-1,5); 对于 FFT和 IFFT的算法和频谱分析的研究 subplot(212); plot(f2,abs(M2); xlabel(f); ylabel(三角波频谱分析后的频谱); axis(-1,1,0,30); fi

36、gure(3) y1=F2T(M1,fs); yy1=y1(1:length(t); subplot(211); plot(t,-yy1); xlabel(t); ylabel(三角函数恢复后的时域函数); axis(-0.05,0.05,-5,5); y2=F2T(M2,fs); yy2=y2(1:length(t); subplot(212); plot(t,abs(yy2); xlabel(t); ylabel(三角波恢复后的时域函数); axis(-10,10,-1,5); 对于 FFT和 IFFT的算法和频谱分析的研究 function M,m,df1,f=T2F(m,ts,df,f

37、s) %-输入参数 %m:信号 %ts:系统时域采样间隔 %df:所需的频率分辨率 %fs:系统采样频率 %-输出(返回)参数 %M:傅里叶变换后的频谱序列 %m: 输入信号参与过傅里叶变换后对应的序列，需要注意的是，该序列与输入信号 m 的区别，其长度是不一样的，输入的 m 长度不一定是 2 的整数次幂，而傅里叶变换要求输入信号长度为 2 的整数次幂，故傅里叶变换前需对输入的 m 信号进行补零，其长度有所增加，故输出参数中的 m 为补零后的输入信号，其长度与输入参数 m 不一样， 但与 M,f 长度是一样的， 并且， 其与时间序列 t 所对应的序列 m(1:length(t)与输入参数中的

38、m 是一致的。 %df1:返回的频率分辨率 %f:与 M 相对应的频率序列 M,m,df1=fftseq(m,ts,df); f = 0:df1:df1*(length(m)-1) -fs/2; %频率向量 M=M/fs; function out=weidaoxu(in1,long) in=in1-1; 对于 FFT和 IFFT的算法和频谱分析的研究 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