对于FFT和IFFT的算法和频谱分析的研究.docx
《对于FFT和IFFT的算法和频谱分析的研究.docx》由会员分享，可在线阅读，更多相关《对于FFT和IFFT的算法和频谱分析的研究.docx（20页珍藏版）》请在冰豆网上搜索。
对于FFT和IFFT的算法和频谱分析的研究
实验总成绩：
——————————————————————————装订线—————————————————————————————————
报告份数：
西安邮电大学
通信与信息工程学院
科研训练论文
专业班级：
学生姓名：
学号（班内序号）:
2013年9月22日
——————————————————————————装订线—————————————————————————————————
对于FFT和IFFT的算法和频谱分析的研究
（ThealgorithmsandspectrumanalysisofFFTandIFFT）
摘要：
快速傅里叶变换（FastFourierTranformation，FFT）是将一个大点数N的DFT分解为若干小点的DFT的组合。
将运算工作量明显降低，从而大大提高离散傅里叶变换（DFT）的计算速度,从而更加适合进行实时运算。
因各个科学技术领域广泛的使用了FFT技术它大大推动了信号处理技术的进步，现已成为数字信号处理强有力的工具，本论文将比较全面的叙述各种快速傅里叶变换算法原理、特点，并完成了基于MATLAB的实现。
最后通过FFT和IFFT的两个应用升余弦滚降和确定函数的频谱分析来分别验证FFT和IFFT的正确性和优越性。
关键词：
MATLABFFTIFFT升余弦滚降函数三角函数三角波函数
Abstract：
FastFourierTransformationismainlyachievedbydevidingthe
DFTwhichisconsistedofNpointsintoalotoflittlesizeofDFT,thustheworkloadcouldbereducedtosomeextent.Itisobviousthatthereduceofworkloadmayimprovethecomputingspeed,soitisfittogettheresultsimmediately.Becauseofitswideapplicationinvariousfieldsofscience,thecomputingspeedisconsiderablyimproved.FFTandIFFTpromotetheadvancementofsignalprocessingespeciallyDigitalSignalProcessing,soithasbecomeapowerfultoolinDSP.Thisarticlewillintroducethetheory,thefeatureofFFTandIFFT,thenwewillgetitachievedbasedonMATLAB.FinallywewillgettwoapplicationsofFFTandIFFTincludingraisedcosineroll-offandthespectrumanalyzingoftwodeterminedfunctiontoverifytheFFTandIFFTcorrectnessandsuperiorityrespectively.
Keywords:
MATLAB,FFT,IFFT,Raisedcosinefunction,Trigonometric, Triangularwavefunction
引言
1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在《计算数学》杂志上发表了“机器计算傅立叶级数的一种算法”的文章，这是一篇关于计算DFT的一种快速有效的计算方法的文章。
它的思路建立在对DFT运算内在规律的认识之上。
这篇文章的发表使DFT的计算量大大减少，并导致了许多计算方法的发现。
这些算法统称为快速傅立叶变换（FastFourierTransform），简称FFT，1984年，法国的杜哈梅尔（P.Dohamel）和霍尔曼（H.Hollmann）提出的分裂基快速算法[2][6]
使运算效率进一步提高。
FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
随着科学的进步，FFT算法的重要意义已经远远超过傅里叶分析本身的应用。
FFT算法之所以快速，其根本原因在于原始变化矩阵的多余行，此特性也适用于傅里叶变换外的其他一些正交变换，例如，快速沃尔什变换、数论变换等等。
在FFT的影响下，人们对于广义的快速正交变换进行了深入研究，使各种快速变换在数字信号处理中占据了重要地位。
因此说FFT对数字信号处理技术的发展起了重大推动作用。
1.FFT的算法
1.1FFT算法的基本思想[1]
设离散的有限长时间序列x（n），
0≤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]
1.3算法的程序实现思想及分析
首先检验待变换的序列的元素个数是否为2的幂次方个，如果不是的话则将其补零使之成为2的整数幂次方个。
然后利用已经编好的位倒序子程序输出位倒序序号，将输入序列不断分组，进行处理后从新分组，直至完成最后的处理即可输出变换后的结果。
需要注意的是，这里fft的变换后结果的元素个数可能与原输入序列的个数不一样，因为如果不是2的幂次方个的话输入序列后面是要补零的。
（主程序为fft_dit和fft_dif）
1.4流程示意图
整个FFT频谱分析与显示过程可用图2所示的流程图示意.
fft按频域抽取算法
fft按时域抽取算法
图2
2.IFFT的频谱变换基本原理
在实验中为了简化算法我们直接利用前面已经编好的fft_dit或fft_dif这两个现成函数来实现，基本原理如下：
将要变换的频谱序列先取共轭然后将其送入前面的函数，将变化后的结果再取共轭即实现IFFT的功能。
主程序为ifft_my
图3
3升余弦滚降
3.1升余弦滚降原理
要实现无码间干扰基带传输时，系统必须满足奈奎斯特准则即：
对于上述公式，我们分3种情况来说明其含义：
（1）Ts<1/2W,其中Ts为系统的输入数据的符号间隔，W为系统的传递函数 X（f）的截止频率。
由于：
因而Z（f）是由频率间隔为1/Ts的X（f）曲线无频率重叠地周期性复制构成。
（2）若Ts=1/2W。
Z（f）仍是由频率间隔为1/Ts的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]OppenheimAV,SchaferRW.DigitalSignalProcessing[M].NY:
Prentice-Hall,1975.
[3]ProakisJG,ManolakisDG.DigitalSignalProcessing-Principles[M].Algorithmsandapplications（2ndEdition）,NY:
Printice-Hall,1995.
[4]崔灵智.Matlab在数字信号处理课程设计中的应用[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）
ifnargin==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）;
forn=1:
N
num（n）=in（weidaoxu（n,N）+1）;
end
forn=1:
k
forc=1:
2^（k-n）
num1=zeros（1,2^n）;
num2=zeros（1,2^n）;
num1=num（（c-1）*2^n+1:
c*2^n）;
form=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）
ifnargin==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）;
forn=1:
k
forc=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））;
form=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
forn=1:
N
num3（n）=num（weidaoxu（n,N）+1）;
end
（3）weidaoxu.m位倒序函数
functionout=niweidaoxu（in,long）
l=log2（long）;
num1=zeros（1,l）;
%out=0;
fori=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来分析升余弦滚降函数的时域和频域的特点
clearall;
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（[-2201.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.250.25-0.10.5]）;
a=0.5;
[f,out_f]=shengyuxiangunjiang（a,N,l）;
subplot（223）;
plot（f,out_f）;
xlabel（'f'）;
ylabel（'a=0.5时的频谱幅度'）;
axis（[-2201.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.250.25-0.10.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（[-2201.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.250.25-0.10.5]）;
a=1;
[f,out_f]=shengyuxiangunjiang（a,N,l）;
subplot（223）;
plot（f,out_f）;
xlabel（'f'）;
ylabel（'a=1时的频谱幅度'）;
axis（[-2201.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.250.25-0.10.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;
form=1:
N/2
out_f（N+1-m）=out_f（m）;
end
（3）rand01.m产生随机的01序列
functions=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主程序
clearall
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）],zeros（1,length（t）/2）];
forn=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）;
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]）;
（2）T2F.m改进后的傅里叶变换[6]
function[M,m,df1,f]=T2F（m,ts,
升级会员 