8点基于DIT的FFT的实现文档格式.docx
《8点基于DIT的FFT的实现文档格式.docx》由会员分享,可在线阅读,更多相关《8点基于DIT的FFT的实现文档格式.docx(11页珍藏版)》请在冰豆网上搜索。
3.6调用系统函数验证10
4心得体会11
参考文献12
摘要
FFT,即为快速傅氏变换,是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。
本文详细介绍了快速傅里叶算法的原理及推详细推导过程,并给出了8点fft的蝶形图及matlab程序代码,并通过调用该函数就算16点的fft。
关键词:
matlab;
fft;
函数
1FFT原理与实现
1.1引言
傅利叶变换是一种将信号从时域变换到频域的变换形式,是声学、语音、电信和信号处理等领域中一种重要的分析工具。
离散傅利叶变换(DFT)是连续傅利叶变换在离散系统中的表示形式,由于DFT的计算量很大,因此在很长一段时间内其应用受到很大的限制。
20世纪60年代由Cooley和Tukey提出了快速傅利叶变换(FFT)算法,它是快速计算DFT的一种高效方法,可以明显地降低运算量,大大地提高DFT的运算速度,运算时间缩短一至两个数量级,从而使DFT在实际应用中得到了广泛的应用。
(k=0,1,..N-1
)
(1)
1.2DFT计算公式
其中x(n)表示输入的离散数字信号序列,WN为旋转因子,X(k)一组N点组成的频率成分的相对幅度。
一般情况下,假设x(n)来自于低通采样,采样频率为fs,那么X(k)表示了从-fs/2率开始,频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。
因为DFT计算得到的一组离散频率幅度只实际上是在频率轴上从成周期变化的,即X(k+N)=X(k)。
因此任意取N个点均可以表示DFT的计算效果,负频率成分比较抽象,难于理解,根据X(k)的周期特性,于是我们又可以认为X(k)表示了从零频率开始,频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。
1.3旋转因子WN的特性
(1)WN的对称性
(2)
(2)WN的周期性
(3)
(3)WN的可约性
(4)
根据以上这些性质,我们可以得到式(5)的一些列有用结果
(5)
假设采样序列点数为N=2^L,L为整数,如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。
首先我们将序列x(n)按照奇偶分为两组如下:
X(2r)=x1(r),X(2r+1)=x2(r),r=0,1..N/2-1(6)
于是根据DFT计算公式
(1)有:
(7)
至此,我们将一个N点的DFT转化为了式(7)的形式,此时k的取值为0到N-1,现在分为两段来讨论,当k为0~N/2-1的时候,因为x1(r),x2(r)为N/2点的序列,因此,此时式(7)可以写为:
k=0...N-1(8)
而当k取值为N/2~N-1时,k用k’+N/2取代,k’取值为0~N/2-1。
对式(7)化简可得:
(9)
综合以上推导我们可以得到如下结论:
一个N点的DFT变换过程可以用两个N/2点的DFT
变换过程来表示,其具体公式如式(10)所示DFT快速算法的迭代公式:
(10)
X1(k)
X2(k)
上式中X'
(k’)为偶数项分支的离散傅立叶变换,X'
'
(k’’)为奇数项分支的离散傅立叶变换。
式(10)的计算过程可以用图1的蝶形算法流图直观地表示出来。
图1时间抽取蝶形运算
在图1中,输入为两个N/2点的DFT输出为一个N点的DFT结果,输入输出点数一致。
运用这种表示方法,8点的DFT可以用图2来表示:
根据公式(10),一个N点的DFT可以由两个N/2点的DFT运算构成,再结合图1的蝶形信号流图可以得到图2的8点DFT的第一次分解。
该分解可以用以下几个步骤来描述:
1.将N点的输入序列按奇偶分为2组分别为N/2点的序列
2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2
3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果
根据式(10)我们可以对图2中的4点DFT进一步分解,得到图3的结果,分解步骤和前面一致。
最后对2点DFT进一步分解得到最终的8点FFT信号计算流图:
从图2到图4的过程中关于旋转系数的变化规律需要说明一下。
看起来似乎向前推一级,
在奇数分组部分的旋转系数因子增量似乎就要变大,其实不是这样。
事实上奇数分组
部分的旋转因子指数每次增量固定为1,只是因为每向前推进一次,该分组序列的数据个
数变少了,为了统一使用以原数据N为基的旋转因子就进行了变换导致的。
每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。
以上边的8点DFT为例,第一次分组N=8,第二次分组N为4,为了统一根据式(4)进行了变换将N变为了8,但指数相应的需要乘以2。
1.4调用8点计算16点
先将16点fft拆成两个8点fft,分别调用自身8点fft函数fft8计算,然后再按以下公式计算出16点的fft
(11)
2程序代码
2.1计算8点FFT代码
functiony=fft8(x)%根据蝶形图计算8点FFT
wn=exp(-j*2*pi/8);
%旋转因子
x1
(1)=x
(1)+x(5);
%计算蝶形图第一级
x1
(2)=x
(1)-x(5);
x1(3)=x(3)+x(7);
x1(4)=x(3)-x(7);
x1(5)=x
(2)+x(6);
x1(6)=x
(2)-x(6);
x1(7)=x(4)+x(8);
x1(8)=x(4)-x(8);
x2
(1)=x1
(1)+x1(3);
%计算蝶形图第二级
x2
(2)=x1
(2)+x1(4)*(wn.^2);
x2(3)=x1
(1)-x1(3);
x2(4)=x1
(2)-x1(4)*(wn.^2);
x2(5)=x1(5)+x1(7);
x2(6)=x1(6)+x1(8)*(wn.^2);
x2(7)=x1(5)-x1(7);
x2(8)=x1(6)-x1(8)*(wn.^2);
y
(1)=x2
(1)+x2(5);
%计算蝶形图输出
y
(2)=x2
(2)+x2(6)*(wn.^1);
y(3)=x2(3)+x2(7)*(wn.^2);
y(4)=x2(4)+x2(8)*(wn.^3);
y(5)=x2
(1)-x2(5);
y(6)=x2
(2)-x2(6)*(wn.^1);
y(7)=x2(3)-x2(7)*(wn.^2);
y(8)=x2(4)-x2(8)*(wn.^3);
2.2计算16点FFT代码
functiony=fft16(x)%16点FFT
wn=exp(-j*2*pi/16);
y1=fft8(x(1:
2:
16));
%计算偶数组的8点fft
y2=fft8(x(2:
%就算奇数组的8点fft
y(1:
8)=y1+y2.*(wn.^[0:
7]);
%计算前八个点
y(9:
16)=y1-y2.*(wn.^[0:
%计算后八个点
3调试过程
3.1编写8点FFT函数
运行matlab,点击file-new-function打开函数编辑窗,输入代码并保存为fft8.m
3.28点FFT函数运行结果
在命令窗口输入fft8(1:
8)后按回车。
3.3调用系统函数验证
在命令窗输入fft(1:
8)后按回车,通过调用系统函数fft计算8点fft,计算结果与上面一致,结果如下:
3.4编写16点FFT函数与运行结果
在8点fft函数的基础上编写16点fft计算函数,步骤如3.3,并在命令窗输入x=[1234567890123456];
fft16(x)后按回车,结果如下
3.5调用系统函数
在命令窗口输入x=[1234567890123456];
fft(x)后按回车,通过调用系统函数fft计算16点fft,计算结果与上面一致,结果如下:
经过验证,所编写的函数正确,可以实现16点的FFT运算。
4心得体会
设计过程中,学习了许多数字信号处理课程中关于数字滤波器的设计的内容,再通过利用参考文献与网络,完成了用Matlab进行数字信号处理课程设计。
通过课程设计,加深了对课堂抽象概念的理解,巩固了课堂上所学的理论知识,并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。
同时掌握编程方法和解决实际问题的技巧。
与其他高级语言的程序设计相比,MATLAB环境下更方便、快捷,节省大量的编程时间,提高编程效率,且参数的修改也十分方便,还可以进一步进行优化设计。
参考文献
[1]邹其洪.《MATLAB教程》.电子工业出版社,2005
[2]吴友宇.《数字信号处理》.东南大学出版社,2008
[3]吴锡龙.《信号与系统》.高等教育出版社,2004
[4]陈怀琛.《MATLAB应用与提高》.西安电子科技大学出版社,2000
[5]丁春利.《DSP技术》.清华大学出版社,2002