MATLAB实验傅里叶分析
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验七 傅里叶变换
一、实验目的
傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。通过傅里叶变换技术可以将时域上的波形分 布变换为频域上的分布,从而获得信号的频谱特性。MA TLAB 提供了专门的函数fft 、ifft 、fft2(即2维快速傅里叶变换)、ifft2以及fftshift 用于实现对信号的傅里叶变换。本次实验的目的就是练习使用fft 、ifft 以及fftshift 函数,对一些简单的信号处理问题能够获取其频谱特性(包括幅频和相频特性)。
二、实验预备知识
1. 离散傅里叶变换(DFT)以及快速傅里叶变换(FFT)简介
设x (t )是给定的时域上的一个波形,则其傅里叶变换为
2()() (1)j ft X f x t e dt π∞
--∞
=⎰
显然X ( f )代表频域上的一种分布(波形),一般来说X ( f )是复数。而傅里叶逆变换定义为:
2()() (2)j ft x t X f e df π∞
-∞
=⎰
因此傅里叶变换将时域上的波形变换为频域上的波形,反之,傅里叶逆变换则将频域上的波形变换为时域上的波形。
由于傅里叶变换的广泛应用,人们自然希望能够使用计算机实现傅里叶变换,这就需要对傅里叶变换(即(1)式)做离散化处理,使之符合电脑计算的特征。另外,当把傅里叶变换应用于实验数据的分析和处理时,由于处理的对象具有离散性,因此也需要对傅里叶变换进行离散化处理。而要想将傅里叶变换离散化,首先要对时域上的波形x (t )进行离散化处理。采用一个时域上的采样脉冲序列:
δ (t -nT ), n = 0, 1, 2, …, N -1;
可以实现上述目的,如图所示。其中N 为采样点数,T 为采样周期;f s = 1/T 是采样频率。注意采样时,采样频率f s 必须大于两倍的信号频率(实际是截止频率),才能避免混迭效应。
接下来对离散后的时域波形()()()()x t x t t nT x nT δ=-=的傅里叶变换()X f 进行离散处理。与上述做法类似,采用频域上的δ脉冲序列:
δ ( f -n/T 0), n = 0, 1, 2, …, N -1;T 0= NT 为总采样时间
可以实现傅里叶变换()X f 的离散化,如下图示。不难看出,离散后的傅里叶变换其频率间隔(频率轴上离散点的间隔,即频域分辨率)
x (t )
δ 脉冲序列
x (t )δ (t -nT )
t
t
t
011
(3)s f f T NT N
∆===
因此要增加分辨率须增加采样点数目N 。频域上每个离散点对应的频率为:
0; 0,1,2,...,-1 (4)s n f n n
f n n N T NT N
=
===显然n = 0的点对应于直流成分。
经过以上离散化处理之后,连续积分的傅里叶变换(1)式转变为如下离散形式:
1
2/0
()(), 0,1,2,..., 1 (5)
N j nk N n k k X f x t e n N π--===-∑其中t k = kT (k =0,1,2,…,N-1)代表采样点时刻。X ( f n )一般是复数,因此离散傅里叶变换(DFT)后变成一个N 点(采样点数)的复数序列。X ( f n )绝对值代表振幅,其幅角代表相位,因此由(5)式可以给出DFT 的振幅频谱和相位频谱。(5)式通常又简写成如下形式:
1
0()(), 0,1,2,..., 1 (6)N nk
N k X n x k W n N -===-∑
其中 2/ j N N W e π-=,x 是采样点数据,它是一个N 个点的向量,DFT 的结果X 是N 个点的复数向量。(5)式或(6)式就是对傅里叶变换进行数值计算的基础。
一般采样点数N 越大,DFT 的结果越接近真实的情况,但是当N 较大时,(6)式的计算量很大,因为使用计算机求解(6)式时,总共要执行N 2次复数乘法和N×(N-1)次复数加法。所以直接用DFT 算法(即(5)式)进行谱分析和信号的实时处理是不切实际的。为了减轻计算的压力,人们提出了一种所谓快速傅里叶变换(FFT )的思想:
取N =2m ,首先将N 个点的采样数据011[,,...,]N x x x x -=分成两个N /2点的序列:
1022[,,...,]N x x x x -= (偶数序列) 2131[,,...,]N x x x x -= (奇数序列)
这样处理的好处是可以把(6)式分解为两个N/2点的DFT ,使计算量降下来。接下来再将N/2点的序列x 1仿照上述做法进一步分裂成2个N/4点的序列x 3和x 4,另一序列x 2亦做如此处理,分裂成2个N/4点的序列x 5和x 6。这样两个N/2点的序列分成了更短的4个N/4点的序列,依次类推,最后的结果是将一个N 点的序列x 裂成了N 个点的单点序列:x 0, x 1, x 2, …, x N-1。这样做可以将DFT 的运算效率提高1-2个数量级,为数字信号处理技术应用于各种信号的实时处理创造了条件,从而推动数字处理技术的发展。由此可见FFT 的思想实质是
f
f
f
()X f
δ (f -n/T 0)
X ( f )
混迭
不断地把长序列的DFT计算分解成若干短序列的DFT,并利用旋转因子(即W N )的周期性和对称性来减少DFT的运算次数。所以FFT就是DFT的快速算法。
有关FFT算法的详细介绍和理论推导参见有关的书籍,这里不做进一步介绍。
2. FFT的MATLAB实现
为了实现快速傅里叶变换,MA TLAB提供了fft、ifft、fft2、ifft2以及fftshift函数,分别用于一维和二维离散傅里叶变换(DFT)及其逆变换。借助这些函数可以完成很多信号处理任务。考虑到信号处理包含的领域很广泛,这里只介绍一维傅里叶变换及其逆变换函数。
(1) fft函数
该函数使用了快速算法来实现时域信号的离散傅里叶变换。常用的格式:
Y= fft (x)
Y= fft (x, m)
Y返回值(复数),返回m点的DFT序列,即(6)式左边的X;
m计算时使用的数据点数(样本数);
x时域信号x(t)在采样点t k处的值,即(6)式右边的x;若实际采样点数目为N(m和N都须是2的幂次),则x为N个元素(即长度N)的向量;若向量x的长度小于m,那么计算时将自动在x序列的后面补0;若x的长度大于m,则x自动截断,使之长度为m。对信号进行频谱分析时,数据样本应有足够的长度,一般FFT程序中所用数据点数(m)最好与原信号含有的数据点数(即输入的样本数N)相同,这样的频谱图具有较高的质量,可减小因补零或截断而产生的影响。
两点说明:
①关于FFT振幅频谱和相位频谱的计算
由于傅里叶变换的结果一般是复数,所以
●对fft的结果取绝对值abs()可以得到振幅,即
Amplitude = abs(Y)
需要注意的是这样得到的幅值实际并非真正的信号振幅,因其值与FFT使用的数据点数N 有关,但不影响分析结果,在IFFT(逆变换)时已经做了处理。要得到真实的振幅值的大小,只要将上述结果除以N/2即可。
●对fft的结果使用函数angle()可以得到相位的结果。但是使用angle函数计算复数的相角时,系统规定一、二象限的角为0~π;三、四象限的角为-π~0。因此若一个角度本来应该从0变到2π,但计算得到的结果却是0~π,再由-π~0,在π处发生跳变,跳变幅度为2π,这就叫相位的卷绕。这种相位的卷绕会使得相频图不连续,呈现锯齿状,为了平滑相频图,通常要再使用unwrap()函数进行相位的解卷绕。因此FFT的相位频谱图应该如下实现Phase = unwrap(angle(Y))
②FFT的振幅频谱具有对称性