matlab课程结题报告.docx
《matlab课程结题报告.docx》由会员分享，可在线阅读，更多相关《matlab课程结题报告.docx（14页珍藏版）》请在冰豆网上搜索。
matlab课程结题报告
《MATLAB与信号处理系统课程设计》
课程性质：
考察
学号：
姓名：
专业：
通信工程
授课教师：
完成日期：
快速傅里叶变换-基2时间抽取FFT算法matlab实现
作者姓名：
摘要：
FFT是一种快速的傅里叶变换。
DFT是信号分析与处理中的一种重要变换。
人们不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数，FFT随之产生了。
MATLAB提供的FFT函数是一个计算DFT的智能程序，能自动选择快速算法进行DFT运算.按基2时间抽取FFT算法，根据傅里叶变换的原理和规律，绘出了算法实现的程序框图，列出了MATLAB环境下软件实现的程序，建立了从算法理论到程序实现的完整概念。
关键字：
离散傅里叶变换;快速傅里叶变换;按时间抽取;MATLAB
DFT是信号分析与处理中的一种重要变换。
但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。
所以快速傅里叶变换FFT在此时就出现了，大大提高了DFT算法的效率，推动了数字信号处理技术的发展。
有限长序列x（n）的N点DFT定义为：
，式中
，其整数次幂简称为旋转因子。
直接进行DFT运算大约需要
次三角函数计算、
次实数乘法计算和
次实数加法计算，且需许多索引和寻址操作［2］。
本文列出了直接DFT的MATLAB程序，这种直接DFT运算概念清楚、编程简单，但占用内存大、运算速度低，在实际工作中并不实用。
基2FFT算法的基本思想是把原始的N点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT，再进行适当的组合，得到原N点序列的DFT，最终达到减少运算次数，提高运算速度的目的。
1.DIT-FFT算法的运算规律
DFT的算法和时间复杂度
对于一个长度为N的离散信号序列x[n]，其DFT变换为
（1）
其中
对任意
完成
（2）式的计算要做N次复数乘法和N-1次复数加法，那么对于
（1）就要做N次复数乘法和N（N-1）次复数加法，其时间复杂度是O（N2），计算工作量和运算时间是非常巨大的。
故要对其进行性质分析，进而做出快速的算法实现DFT，以便计算机能够快速实现，从而在实际信号处理工作中得到广泛的应用。
按时间抽取的基2FFT算法，先是将N点输入序列x（n）在时域按奇偶次序分解成2个N/2点序列x1（n）和x2（n），再分别进行DFT运算，求出与之对应的X1（k）和X2（k），然后利用下图所示的运算流程进行蝶形运算，得到原N点序列的DFT。
用N/2点X1（k）和X2（k）表示序列x（n）的N点DFTX（k）
注意：
这里的k的取值范围为0,1，…N-1
蝶形运算符号
DIT蝶形运算流图符号
为了编写DIT-FFT算法的运算程序，首先要分析其运算规律，总结编程思想并绘出程序框图。
由下图可知，DIT-FFT算法的运算过程很有规律。
8点DIT-FFT运算流程
2．编程思想
实现FFT运算的核心是蝶形运算，找出蝶形运算的规律是编程的基础。
蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。
对
点的FFT共有M级运算，用L表示从左到右的运算级数（L=1，2，…，M）。
第L级共有
个不同指数的旋转因子，用R表示这些不同指数旋转因子从上到下的顺序（R=0，1，…，B-1）。
第R个旋转因子的指数
，旋转因子指数为P的第一个蝶的第一节点标号k从R开始，由于本级中旋转因子指数相同的蝶共有
个，且这些蝶的相邻间距为
，故旋转因子指数为P的最后一个蝶的第一节点标号k为：
，本级中各蝶的第二个节点与第一个节点都相距B点。
应用原位计算，蝶形运算可表示成如下形式：
图3DIT-FFT运算程序框图
总结上述运算规律，可采用如下运算方法进行DIT-FFT运算。
首先读入数据，根据数据长度确定运算级数M，运算总点数
，不足补0处理。
然后对读入数据进行数据倒序操作。
数据倒序后从第1级开始逐级进行，共进行M级运算。
在进行第L级运算时，先算出该级不同旋转因子的个数
（也是该级中各个蝶形运算两输入数据的间距），再从R=0开始按序计算，直到R=B-1结束。
每个R对应的旋转因子指数
，旋转因子指数相同的蝶从上往下依次逐个运算，各个蝶的第一节点标号k都是从R开始，以
为步长，到
（可简取极值N-2）结束。
考虑到蝶形运算有两个输出，且都要用到本级的两个输入数据，故第一个输出计算完毕后，输出数据不能立即存入输入地址，要等到第二个输出计算调用输入数据完毕后才能覆盖。
这样数据倒序后的运算可用三重循环程序实现。
整个运算流程如图所示。
3.程序源代码
%基2DIT-FFT运算的MATLAB程序
clc;closeall;clear;formatcompact;
%输入数据并计算常量
xn=[0,1,2,3,4,5,6,7];%可取任意序列
M=nextpow2（length（xn））,N=2^M,
form=0:
N/2-1;%旋转因子指数范围
WN（m+1）=exp（-j*2*pi/N）^m;%计算旋转因子
end
A=[xn,zeros（1,N-length（xn））];%数据输入
disp（'输入到各存储单元的数据:
'）,disp（A）;
%数据倒序操作
J=0;%给倒序数赋初值
forI=0:
N-1;%按序交换数据和算倒序数
ifIT=A（I+1）;A（I+1）=A（J+1）;A（J+1）=T;
end
%算下一个倒序数
K=N/2;
whileJ>=K;
J=J-K;K=K/2;
end
J=J+K;
end
disp（'倒序后各存储单元的数据:
'）,disp（A）;
%分级按序依次进行蝶形运算
forL=1:
M;%分级计算
disp（'运算级次:
'）,disp（L）;
B=2^（L-1）;
forR=0:
B-1;%各级按序蝶算
P=2^（M-L）*R;
forK=R:
2^L:
N-2;%每序依次计算
T=A（K+1）+A（K+B+1）*WN（P+1）;
A（K+B+1）=A（K+1）-A（K+B+1）*WN（P+1）;
A（K+1）=T;
end
end
disp（'本级运算后各存储单元的数据:
'）,disp（A）;
end
disp（'输出各存储单元的数据:
'）,Xk=A,
disp（'调用fft函数运算的结果:
'）,fftxn=fft（xn,N）,
3.实验及结果分析（核心部分）
题目：
快速傅里叶变换-基2时间抽取FFT算法matlab实现
要求：
已知输入信号x（t）=0.6sin（200πt）+sin（400πt）+0.3sin（800πt）。
1、实现基2时间抽取的FFT算法的1024点，4096点变换，并与matlab自带函数进行比较，包括运算时间和精度。
解：
快速傅里叶变换程序代码：
functiony=myfft（xr,n）
p=0:
n-1;
nu=log2（n）;
p1=p;
b=zeros（1,n）;
fort=1:
nu;
p2=floor（p1/2）;
b=b*2+（p1-2*p2）;
p1=p2;
end;
yr（p+1）=xr（b+1）;
xr=yr;
t=0:
n/2-1;
forv=0:
n/2-1;
w=exp（-2*i*pi*t/n）;
end;
yr（p+1）=xr（b+1）;
xr=yr;
t=0:
n/2-1;
forv=0:
n/2-1;
w=exp（-2*i*pi*t/n）;
end;
form=1:
nu;
h=2^（m-1）;
k=1;
while（kfort=1:
h;
y=bitshift（k-1,nu-m,nu）+1;
xch（k）=xr（k）+w（y）*xr（k+h）;
k=k+1;
end;
fort=1:
h;
y=bitshift（k-1-h,nu-m,nu）+1;
xch（k）=xr（k-h）-xr（k）*w（y）;
k=k+1;
end;
end;
xr=xch;
end;
y=xr
N=1024的傅里叶变换
fs=100;N=1024;
n=0:
N-1;t=n/fs;
x=0.6*sin（200*pi*t）+sin（400*pi*t）+0.3*sin（800*pi*t）;
y=myfft（x,N）;
mag=abs（y）;
f=n*fs/N;
subplot（2,2,1）,plot（f,mag）;
xlabel（'ƵÂÊ/Hz'）;
ylabel（'Õñ·ù'）;title（'N=128'）;gridon;
Subplot（2,2,2）;plot（f,angle（y））;
xlabel（'ƵÂÊ/Hz'）;
ylabel（'Ïàλ'）;title（'N=128'）;gridon;
N=1024的傅里叶自带函数
fs=100;
N=1024;
n=0:
N-1;t=n/fs;
x=0.6*sin（200*pi*t）+sin（400*pi*t）+0.3*sin（800*pi*t）;
y=fft（x,N）;
mag=abs（y）;
f=（0:
length（y）-1）'*fs/length（y）;
subplot（221）;plot（f,mag）;
xlabel（'Frequence（Hz）'）;ylabel（'Magnitude'）;
title（'N=1024'）;grid;
subplot（222）;plot（f,angle（y））;
xlabel（'Frequence（Hz）'）;ylabel（'Magnitude'）;
title（'N=1024'）;grid;
N=4096的傅里叶变换
fs=100;N=4096;
n=0:
N-1;t=n/fs;
x=0.6*sin（200*pi*t）+sin（400*pi*t）+0.3*sin（800*pi*t）;
y=myfft（x,N）;
mag=abs（y）;
f=n*fs/N;
subplot（2,2,1）,plot（f,mag）;
xlabel（'ƵÂÊ/Hz'）;
ylabel（'Õñ·ù'）;title（'N=4096'）;gridon;
Subplot（2,2,2）;plot（f,angle（y））;
xlabel（'ƵÂÊ/Hz'）;
ylabel（'Ïàλ'）;title（'N=4096'）;gridon;
N=4096的傅里叶变换自带函数
fs=100;
N=4096;
n=0:
N-1;t=n/fs;
x=0.6*sin（200*pi*t）+sin（400*pi*t）+0.3*sin（800*pi*t）;
y=fft（x,N）;
mag=abs（y）;
f=（0:
length（y）-1）'*fs/length（y）;
subplot（221）;plot（f,mag）;
xlabel（'Frequence（Hz）'）;ylabel（'Magnitude'）;
title（'N=4096'）;grid;
subplot（222）;plot（f,angle（y））;
xlabel（'Frequence（Hz）'）;ylabel（'Magnitude'）;
title（'N=4096'）;grid;
.与matlab自带函数进行比较结果（运算时间和精度）:
myfft（）函数精度与matlab自带的函数相比所用的时间比matlab相对来说少一些，精确度要高一些，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。
1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。
如果要提高频率分辨力，则必须增加采样点数，也即采样时间。
频率分辨率和采样时间是倒数关系。
所以采样点数越多，精度越高，越能反映快速傅里叶变换的效率高。
课程设计总结：
FFT算法是在实践应用中广泛的成熟的算法，MATLAB提供的FFT函数可以自动选择快速算法进行DFT计算。
在理论学习的基础上，通过本次的试验，加深了对快速傅里叶变换的理解，熟悉了FFT算法及其程序编写，了解了应用FFT进行信号傅里叶频谱过程中可能出现的问题，一边实际中正确应用FFT.这次课程设计是我对MATLAB语言的理解与应用能力得到了较大的提高，也是我认识到只要深入学习，我们有很大的提升空间。
在课程设计的过程中，我们会遇到很多的问题，如：
编写程序中出现的语法问题等。
我们可以通过查阅书籍，询问老师，与同学交流，在网上查找等方式解决这些问题，是我们解决问题的能力有所提升。
参考文献
[1]余成波，陶红艳，《数字信号处理及MATLAB实现》，北京：
清华大学出版社，2008
[2]楼顺天，李博菡，《基于MATLAB的系统分析与设计》，西安：
西安电子科技大学出版社，1998
[3]曹弋，赵阳，《MATLAB实用教程》，北京：
电子工业出版社，2007
[4]高西全，丁玉美，《数字信号处理》第三版，西安电子科技大学出版社，2000
升级会员 