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

FFT 和IFFT 算法及频谱分析 科研训练论文

西安邮电大学

通信与信息工程学院

科研训练论文

专业班级： 学生姓名： 学号(班内序号):

2013 年 9 月 22日

——————————————————————————装

订

线

——————————————————————————————————————————————————

报告份数：

实验总成绩：

对于FFT和IFFT的算法和频谱分析的研究（The algorithms and spectrum analysis of FFT

and IFFT）

摘要：

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

关键词：

MATLAB FFT IFFT 升余弦滚降函数三角函数三角波函数

Abstract：

Fast Fourier Transformation is mainly achieved by deviding the DFT which is consisted of N points into a lot of little size of DFT,thus the

workload could be reduced to some extent.It is obvious that the reduce of

workload may improve the computing speed,so it is fit to get the results

immediately.Because of its wide application in various fields of science , the

computing speed is considerably improved.FFT and IFFT promote the advancement of signal processing especially Digital Signal Processing,so it has become a powerful tool in DSP.This article will introduce the theory,the feature of FFT and IFFT,then we will get it achieved based on MATLAB.Finally we will get two applications of FFT and IFFT including raised cosine roll-off and the spectrum analyzing of two determined function to verify the FFT and IFFT

correctness and superiority respectively.

Key words:

MATLAB ,FFT ,IFFT,Raised cosine function,Trigonometric,

Triangular wave function

引言

1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在《计算数学》杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。它的思路建立在对DFT运算内在规律的认识之上。这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。这些算法统称为快速傅立叶变换（Fast Fourier Transform），简称FFT，1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法[2][6]

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

1.FFT的算法

1.1 FFT算法的基本思想[1]

设离散的有限长时间序列x(n)，

0≤n≤N-1，则其离散傅立叶变换为：

这样，矩阵W 中有许多相同的元素，从而可

以简化DFT的运算过程．FFT 算法有许多形式，

笔者只讨论最基本的时间抽取基-2FFT 算法．

1.2 算法分析

一个N 点长序列，直接用DFT 方法需要复数乘法N2次；复数加法N(N-1)次。而由

图2 可知，采用FFT 则只需要复数乘法次；复数加法次。当

时，这样，运算速度提高了1-2 个数量级．图1 为FFT 算法和直接DFT 算

法所需运算量与计算点数N 的关系曲线．显然，N 越大时，优越性越明显．但当

N 相当大时，利用单机串行进行FFT 运算同样满足不了实时系统的需要．[ 1 ]

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(6) x(7)

X(5)

X(4)

X(1)

X(2)

X(6)

X(3)

X(7) fft按频域抽取算法

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

待变换序列取共轭

快速傅里叶变换

取共轭

快速逆傅里叶变换的结

果

图3

3 升余弦滚降

3.1升余弦滚降原理

要实现无码间干扰基带传输时，系统必须满足奈奎斯特准则即：

()m m m X f Ts Ts =∞

=-∞

+=∑ 对于上述公式，我们分3种情况来说明其含义：

（1）Ts<1/2W,其中Ts 为系统的输入数据的符号间隔，W 为系统的传递函数 X （f ）的截止频率。由于：

()()

m m m Z f X f Ts

=∞

=-∞

=

+∑

因而Z （f ）是由频率间隔为1/Ts 的X （f ）曲线无频率重叠地周期性复制构成。 （2）若Ts=1/2W 。Z （f ）仍是由频率间隔为1/Ts 的X （f ）曲线无频率重叠地周期性复制构成，在此情况下，仅有一个情况可满足无码间干扰传输的条件，即当

,0,()Ts f W

X f ?≤?

?=其他

此基带传输系统的传递函数是理想低通，其频带宽度为W，则该系统无码间干扰传输的最小Ts=1/2W,即无码间干扰传输的最大符号速率Rs=1/Ts=2W,称此传输速率为奈奎斯特速率。

在此理想情况下，虽然系统的频带利用率达到极限，但是此时x(t)是sinc函数，她是非因果的，是物理不可实现的。并且，此x(t)冲击脉冲形状收敛到0的速度极慢，若在收端低通滤波器输出端的采样时科存在定时误差，则在实际采样时刻的采样值会存在码间干扰

（3）对于Ts>1/2W情况，Z（f）由频率间隔为1/Ts的X（f）曲线无频率重叠地周期性复制并相加构成的，它还是周期性频谱。在这种情况下，有一特定频谱可满足无码间干扰传输的条件，它就是已获广泛应用的升余弦谱。升余弦滤波器的传递函数表示式为：

称α为滚降因子，取值为0≦α≦1。

在α=0时，滤波器的带宽W为1/(2Ts),称为奈奎斯特带宽；α=0.5时，滤波器的截止频率W=（1+α）/（2Ts）=0.75Rs; α=1时，滤波器的截止频率W=Rs

[2][3][4][5]

3.2 升余弦滚降的实现及其与快速傅里叶逆变换的关系

升余弦滚降函数是在基带无码间干扰传输中经常用到的频域函数。其主要特性是升余弦滚降函数经过频域平移叠加后能够成一个在各个频域幅度恒定的频域函数。我们在实现升余弦滚降时可以首先在频域实现α=0，α=0.5，α=0.75和α=1四种升余弦滚降函数，然后通过自己编写的ifft_my进行傅里叶逆变换，并将变换后的时域结果反映在图上，分别对比α不同是对应的时域上的不同的时域特性。

3.3 变换前的时域特性和变换后的频域特性

图4

4 利用fft和ifft进行具体函数频谱分析的实例

4.1 三角函数的频谱分析及其信号的恢复

在实际信号的分析中，三角函数是非常常见和基本的信号，在这里我们就对三角函数进行分析。在分析中我们会碰到要分析的函数的采样点数不是2的整数次幂个，我们要对它进行补零处理。

4.2 三角波函数的分析

在分析中因为可能会用到自己编写的T2F函数和F2T函数，但是这两个函数仍然是建立在自己编写的fft_dit和ifft_my的基础上的。只是对输入变量进行了更加完善的处理，一个是补零，另一个是对频域进行了搬移，将原来pi~2pi的部分

搬到-pi~0，因为大家在看频谱时比较习惯频谱是关于零对称的。

4.3 分析结果

5 结论

5.1 关于fft和ifft函数的分析

在fft和ifft的实现过程中，的确能降低运算的次数，但是也正好印证了一个很著名的理论，“时间换空间，空间换时间”。在fft和ifft的算法中，我们降低运算的次数是以占用更多的空间换来的。每次将要变换的的序列进行分组，然后对每个组进行处理，虽然降低了运算次数，但是也增加了运算空间的占用。5.2 关于升余弦滚降函数的结论

通过图形我们可以观察出，对于不同α的函数，时域主要部分占用的宽度是一样的。但是随着α的增大，时域的起伏是越来越小的。因此我们可以认为，α越大，时域起伏越小，因此在非线性系统中传输时的失真越小。但是与此同时带来的缺点是占用的α越大，频域占用的带宽是越大的，越利于在限带系统中传输，频率的利用率也就越低。

5.3 关于三角函数和三角波函数的频谱分析结论

（1）根据理论分析，三角函数的频谱因该是几个冲激函数的组合，三角波函数的频谱应该是Sa函数平方的形式。而在实际分析中我们通过观察频谱图可以看出分析的结果基本符合理论分析结果。这也侧面印证了我们自己编写的fft_dit和ifft_my的正确性和有效性。

（2）在通过频谱恢复时域时，我们观察到恢复后的三角函数的时域函数有所失真，三角波函数基本无失真。原因可能是在T2F和F2T

函数中可能因为补零改变了元素个数，因此引起了恢复后的函数时域出现失真。

附录1：

[1]蒋冬初,何飞.FFT算法的并行处理研究.湖南城市学院学报(自然科学版) .2005-06-30

[2]Oppenheim A V, Schafer R W.Digit al Signal Proces sing [M] . N Y: Prent ice-Hall , 1975.

[3]Proak is J G, Manolakis D G. Digital Signal Proces sing-Principles[M] . Algorithms and applicat ions ( 2nd Edit ion) , N Y:Print ice-Hall, 1995.

[4]崔灵智. M atlab 在数字信号处理课程设计中的应用[ J]. 山东水利职业学院刊, 2008. 3: 11- 12.

[5]樊昌信曹丽娜通信原理（第六版）国防工业出版社 2006.8

[6]卢光跃黄庆东包志强数字信号处理处理及应用人民邮电出版社 2012.6

[7]王立宁乐光新詹菲 MATLAB与通信仿真人民邮电出版社 2000.01.01

小组成员：

组长：08谷群（03101045）

组员：31王敏（03101069）

03陈若寰（03101040）

04边钰涵（03101041）

任务分配：

资料查阅：王敏

程序设计：谷群

文本编辑：边钰涵

技术校对：陈若寰

附录：

程序代码

1 .自己实现的fft

(1)fft_dit.m 按时域抽取快速傅里叶变换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);

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,2^n);

num2=zeros(1,2^n);

num1=num((c-1)*2^n+1:c*2^n);

for m=1:2^(n-1)

num2(m)=num1(m)+num1(m+2^(n-1))*exp(-i*(m-1)*2*pi/2^n);

num2(m+2^(n-1))=num1(m)-num1(m+2^(n-1))*exp(-i*(m-1)*2*pi/2^n); end

num((c-1)*2^n+1:c*2^n)=num2;

end

end

num3=num;

(2)fft_dif.m 按频域抽取快速傅里叶变换

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);

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)))*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

（3）weidaoxu.m 位倒序函数

function out=niweidaoxu(in,long)

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

（4）ifft_my.m 快速傅里叶逆变换

function [out]=ifft_my(in)

num1=conj(in);

num2=fft_dit(num1);

out=conj(num2)/length(in);

2.升余弦滚降函数的时域频域分析

(1)主函数

%通过FFT和IFFT来分析升余弦滚降函数的时域和频域的特点clear all;

a=0;

N=128;

l=2;

figure(1)

[f,out_f]=shengyuxiangunjiang(a,N,l);

subplot(221);

plot(f,out_f);

xlabel('f');

ylabel('a=0时的频谱幅度');

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]=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)]; subplot(224);

plot(f,out_t_);

xlabel('t');

ylabel('a=0.5时的时域幅度');

axis([-0.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_);

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(f,out_t_);

xlabel('t');

ylabel('a=1时的时域幅度');

axis([-0.25 0.25 -0.1 0.5]);

（2）shengyuxiangunjiang.m 产生升余弦滚降频谱的函数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;

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

(3)rand01.m 产生随机的01序列

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).

3 确定函数的频谱分析

（1） pinpufenxi.m 主程序

clear all

dt=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);

xlabel('t');

ylabel('三角函数频谱分析前的时域');

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(t):5-5*4/length(t)],zero s(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]);

subplot(212);

plot(f2,abs(M2));

xlabel('f');

ylabel('三角波频谱分析后的频谱');

axis([-1,1,0,30]);

figure(3)

[y1]=F2T(M1,fs);