实验三 FFT算法的应用

实验三 FFT算法的应用

一、实验目的

1. 2. 3.

通过实验加深对快速傅立叶变换（FFT）基本原理的理解。 掌握FFT的用于信号的谱分析； 掌握利用FFT计算卷积。

二、实验仪器设备

PC机 MATLAB软件

三、实验原理

离散傅里叶变换（DFT）和卷积是信号处理中两个最基本也是最常用的运算，它们涉及到信号、系统的分析与综合这一广泛的信号处理领域。实际上卷积与DFT之间有着互通的联系：卷积可化为DFT来实现，其它的许多算法，如相关、滤波和谱估计等都可化为DFT来实现，DFT也可化为卷积来实现。 1.MATLAB中DFT的FFT实现

对N点序列x(n)，其DFT变换对定义为：

N?1?nk?X(k)??x(n)WN?n?0?N?11?nk?x(n)?X(k)WN??Nk?0?k?0,1,...,N?1,WN?n?0,1,...,N?1e?j2?/N

显然，求出N点X(k)需要N次复数乘法，N(N-1)次复数加法。众所周知，实现一次复

数乘需要四次实数乘和两次实数加，实现一次复数加则需要两次实数加。当N很大时，其计算量是相当可观的。例如，若N=1024，则需要1,048,576次复数乘法，即4,194,304次实数乘法。所需时间过长，难于“实时”实现。若是处理二维或三维图像，所需的计算量更是大得惊人。

FFT是DFT的一种快速算法。FFT算法将长序列的DFT分解为短序列的DFT。N点的DFT先分解为2个N/2点的DFT，每个N/2点的DFT又分解为2个N/4点的DFT。按照此规律，最小变换的点数即所谓的“基数（radix）”。因此，基数为2的FFT算法的最小变换（或称蝶形）是2点DFT。一般地，对N点FFT，对应于N个输入样值，有N个频域样值与之对应。

时析型FFT是将序列逐次奇偶对分，奇数号排成一子序列，偶数号排成一子序列，各子序列的长度为N/2，用同样的方法再将两个子序列又分成长度为N/4的4个子序列，依次这样分下去，最后的子序列长度为1，即一个数，而一个数的DFT就是本身。然后找出子序列的DFT与原序列DFT的关系，就可以合成原序列的DFT。

1)FFT的实现：计算矢量或矩阵的FFT：格式y＝fft（x） y＝fft（x，n）。其中y＝fft（x）是利用FFT算法计算矢量x的离散傅立叶变换，若x为矩阵，则y是对矩阵的每一列向量进行FFT；当x的长度是2 的整数次幂时，则fft函数采用基2时析型FFT，否则采用混合基算法。而y＝fft（x，n）执行n点FFT，当x的长度小于n时，在x尾部补零至n

2

点；当x的长度大于n时，则截断x使其长度为n。

2）IFFT实现：格式y＝ifft（x） y＝ifft（x，n）。用于计算矢量或矩阵的IFFT。 例1：MATLAB仿真程序

%*****************************************************************% %mode--信号的种类。1--正弦波;2--方波;3--锯齿波 %Nfft--FFT点数

%*****************************************************************% function [x]=ffts(mode,Nfft) n=0:Nfft-1;

if mode==1 x=sin(2*pi*n/Nfft);end if mode==2 x=square(2*pi*n/Nfft);end if mode==3 x=sawtooth(2*pi*n/Nfft);end set(gcf,'menubar',menubar); subplot(2,1,1);stem(n,x);

axis([0 Nfft-1 1.1*min(x) 1.1*max(x)]); xlabel('Points-->');ylabel('X(n)'); y=abs(fft(x,Nfft)); subplot(2,1,2);stem(n,y);

axis([0 Nfft-1 1.1*min(y) 1.1*max(y)]); xlabel('frequency--->');ylabel('|X(k)|');

%*****************************************************************%

假设需观察正弦波信号的频谱，对正弦波信号作64点的FFT，则只需在MATLAB的命令窗口下键入：

x=ffts(1,64)

被噪声污染的信号，很难看出它包含的频率分量。通常通过FFT分析可得到其频率成分。 例2：一个50Hz和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，采样频率为1khz。现分析其信号频率成分。

t=0:0.001:0.6;%采样频率为1kHz，时间间隔为0.001 x=sin(2*pi*50*t)+sin(2*pi*120*t);

y=x+1.5*randn(1,length(t));%受污染的正弦信号 subplot(2,1,1); plot(t,y);

xlabel('t');ylabel('y(t)'); Y=fft(y,512);

P=Y.*conj(Y)/512;%计算功率谱密度

f=1000*(0:255)/512;%取结果的一半 subplot(2,1,2); plot(f,P(1:256));

xlabel('f');ylabel('P(f)'); 2.利用FFT快速计算卷积

工程实际需要解决的卷积：y(n)?x(n)?h(n)，但其计算量很大。

而圆卷积为：y(n)?x(n)?h(n)，便于采用FFT算法, 故计算速度快。若将线卷积的两个序列用增补零的方法将长度取为一致，此时两序列的离散线卷积和圆周卷积结果是相等的，这样就则可以通过圆卷积来快速计算线卷积。

利用FFT 计算卷积的步骤

（1）设两序列原长度分别为：N和M，将长度增加到L（2）用FFT法求加长序列的DFT频谱； （3）计算两序列DFT频谱的乘积；

（4）用IFFT求DFT频谱乘积的逆变换，便得两序列的离散线卷积。 例3: 求两矩形序列的卷积：

?1,x(n)???0,0?n?4其它?N?M?1（L

为2的整数次幂）；

和

?1/2,h(n)???0,0?n?5 其它x=[1 1 1 1 1 ];

h=[1/2 1/2 1/2 1/2 1/2 1/2]; L=length(x)+length(h)-1;

L=ceil(log2(L));%B = ceil(A) rounds the elements of A to the nearest integers greater than or

%equal to A. For complex A, the imaginary and real parts are rounded independently.

L=2^L;

subplot(311);stem(x,'.'); xlabel('n');ylabel('x(n)'); axis([1,L,0,1.5]); subplot(312);stem(h,'.'); xlabel('n');ylabel('h(n)'); axis([1,L,0,1.5]); X=fft(x,L); H=fft(h,L); Y=X.*H;

y=ifft(Y);

subplot(313);stem(y,'.'); xlabel('n');ylabel('h(n)') ;

上述快速卷积法适合两序列长度相当的情况，当一个序列为长序列，一个序列是短序列时，常采用分段快速卷积法。即N1??N2，通常不允许等x(n)全部采集齐后再进行卷积，否则使输出相对于输入有较长的延时，另外，若N1?N2?1太大，h(n)要补上太多的零点，很不经济，且FFT的计算时间也要很长。为此，采用分段卷积的方法，即把x(n)分成长度与h(n)相仿的一段段，分别求出每段卷积的结果，然后用相应的方式把它们结合起来，便是总的输出。分段卷积方法主要有两种，即重叠相加法和重叠保留法。具体内容请参考教材相应部分，本实验这部分不作重点要求。

四、实验内容

1.参考例1实现余弦信号、三角波信号的FFT。

2. 一个40Hz和100Hz正弦信号构成的信号，受零均值随机噪声的干扰，采样频率为500Hz。现分析其信号频率成分。 3.利用快速卷积法，计算x(n)???1,?0,0?n?15其它 和

?0.8n,h(n)???0,0?n?15其它的卷积。

五、实验步骤

1．启动MATLAB

2．建立新的m-file文件，输入MATLAB代码。

3．按F5或文件栏的运行，分析结果。结果图像贴入实验报告。

六、实验报告要求

在实验报告中应包含以下内容：

（1）实验目的（2）实验仪器设备，（3）实验原理，（4）实验内容和操作步骤，（5）结果和分析。

七、思考题

1. 2. 3.

DFT和信号频谱之间有什么关系？

N同时取8、16、32时，线性卷积和循环卷积的结果有何不同，为什么？ 分析直接计算线性卷积和利用FFT计算线性卷积的时间。

本文来源：https://www.bwwdw.com/article/tva8.html