信号与系统实验报告 实验4非周期信号的傅里叶变换实验.docx

信号与系统实验报告 实验4非周期信号的傅里叶变换实验.docx

《信号与系统实验报告 实验4非周期信号的傅里叶变换实验.docx》由会员分享，可在线阅读，更多相关《信号与系统实验报告 实验4非周期信号的傅里叶变换实验.docx（11页珍藏版）》请在冰点文库上搜索。

信号与系统实验报告实验4非周期信号的傅里叶变换实验

信号与系统

实验报告

实验四非周期信号的傅里叶变换

实验四 非周期信号的傅里叶变换

一、实验目的

傅里叶变换是通信系统、图像处理、数字信号处理以及物理学等领域内的一种重要的数学分析工具。

通过傅里叶变换技术可以将时域上的波形分布变换为频域上的分布，从而获得信号的频谱特性。

MATLAB提供了专门的函数fft、ifft、fft2（即2维快速傅里叶变换）、ifft2以及fftshift用于实现对信号的傅里叶变换。

本次实验的目的就是练习使用fft、ifft以及fftshift函数，对一些简单的信号处理问题能够获取其频谱特性（包括幅频和相频特性）。

二、实验预备知识

1.离散傅里叶变换（DFT）以及快速傅里叶变换（FFT）简介

设x（t）是给定的时域上的一个波形，则其傅里叶变换为

显然X（f）代表频域上的一种分布（波形），一般来说X（f）是复数。

而傅里叶逆变换定义为：

因此傅里叶变换将时域上的波形变换为频域上的波形，反之，傅里叶逆变换则将频域上的波形变换为时域上的波形。

由于傅里叶变换的广泛应用，人们自然希望能够使用计算机实现傅里叶变换，这就需要对傅里叶变换（即

（1）式）做离散化处理，使之符合电脑计算的特征。

另外，当把傅里叶变换应用于实验数据的分析和处理时，由于处理的对象具有离散性，因此也需要对傅里叶变换进行离散化处理。

而要想将傅里叶变换离散化，首先要对时域上的波形x（t）进行离散化处理。

采用一个时域上的采样脉冲序列:

（t－nT）,n=0,1,2,…,N－1；

可以实现上述目的，如图所示。

其中N为采样点数，T为采样周期；fs=1/T是采样频率。

注意采样时，采样频率fs必须大于两倍的信号频率（实际是截止频率），才能避免混迭效应。

接下来对离散后的时域波形

的傅里叶变换

进行离散处理。

与上述做法类似，采用频域上的脉冲序列：

（f－n/T0）,n=0,1,2,…,N－1；T0=NT为总采样时间

可以实现傅里叶变换

的离散化，如下图示。

不难看出，离散后的傅里叶变换其频率间隔（频率轴上离散点的间隔，即频域分辨率）

因此要增加分辨率须增加采样点数目N。

频域上每个离散点对应的频率为:

显然n=0的点对应于直流成分。

经过以上离散化处理之后，连续积分的傅里叶变换

（1）式转变为如下离散形式：

其中tk=kT（k=0,1,2,…,N-1）代表采样点时刻。

X（fn）一般是复数，因此离散傅里叶变换（DFT）后变成一个N点（采样点数）的复数序列。

X（fn）绝对值代表振幅，其幅角代表相位，因此由（5）式可以给出DFT的振幅频谱和相位频谱。

（5）式通常又简写成如下形式：

其中

，x是采样点数据，它是一个N个点的向量，DFT的结果X是N个点的复数向量。

（5）式或（6）式就是对傅里叶变换进行数值计算的基础。

一般采样点数N越大，DFT的结果越接近真实的情况，但是当N较大时，（6）式的计算量很大，因为使用计算机求解（6）式时，总共要执行N2次复数乘法和N×（N-1）次复数加法。

所以直接用DFT算法（即（5）式）进行谱分析和信号的实时处理是不切实际的。

为了减轻计算的压力，人们提出了一种所谓快速傅里叶变换（FFT）的思想：

取N=2m，首先将N个点的采样数据

分成两个N/2点的序列：

（偶数序列）

（奇数序列）

这样处理的好处是可以把（6）式分解为两个N/2点的DFT，使计算量降下来。

接下来再将N/2点的序列x1仿照上述做法进一步分裂成2个N/4点的序列x3和x4，另一序列x2亦做如此处理，分裂成2个N/4点的序列x5和x6。

这样两个N/2点的序列分成了更短的4个N/4点的序列，依次类推，最后的结果是将一个N点的序列x裂成了N个点的单点序列：

x0,x1,x2,…,xN-1。

这样做可以将DFT的运算效率提高1－2个数量级，为数字信号处理技术应用于各种信号的实时处理创造了条件，从而推动数字处理技术的发展。

由此可见FFT的思想实质是不断地把长序列的DFT计算分解成若干短序列的DFT，并利用旋转因子（即WN）的周期性和对称性来减少DFT的运算次数。

所以FFT就是DFT的快速算法。

有关FFT算法的详细介绍和理论推导参见有关的书籍，这里不做进一步介绍。

2.FFT的MATLAB实现

为了实现快速傅里叶变换，MATLAB提供了fft、ifft、fft2、ifft2以及fftshift函数，分别用于一维和二维离散傅里叶变换（DFT）及其逆变换。

借助这些函数可以完成很多信号处理任务。

考虑到信号处理包含的领域很广泛，这里只介绍一维傅里叶变换及其逆变换函数。

（1）fft函数

该函数使用了快速算法来实现时域信号的离散傅里叶变换。

常用的格式：

Y=fft（x）

Y=fft（x,m）

Ｙ　返回值（复数），返回m点的DFT序列，即（6）式左边的X；

m　计算时使用的数据点数（样本数）；

x　时域信号x（t）在采样点tk处的值，即（6）式右边的x；若实际采样点数目为Ｎ（m和Ｎ都须是２的幂次），则ｘ为Ｎ个元素（即长度Ｎ）的向量；若向量ｘ的长度小于m，那么计算时将自动在ｘ序列的后面补０；若ｘ的长度大于m，则ｘ自动截断，使之长度为m。

对信号进行频谱分析时，数据样本应有足够的长度，一般FFT程序中所用数据点数（m）最好与原信号含有的数据点数（即输入的样本数Ｎ）相同，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。

两点说明：

①关于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的振幅频谱具有对称性

如下图所示。

因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内（共N/2+1个频率点）的幅频特性。

（2）fftshift函数

其作用是将零频点移到频谱的中间（即Nyquist频率处），使用格式：

Y=fftshift（X）

X是向量，该命令将零频点移动到频谱X的中间，并交换频谱X的左右两半。

将零频点放到频谱的中间对于观察傅立叶变换是有用的。

例１：

对时域信号

进行频谱分析。

fs=100;%采样频率>2倍的信号频率

N=256;%采样点数目（=2的幂次）

n=0:

N-1;%构造采样点序列

t=n/fs;%得到采样时间序列，t=nT=n/fs

x=0.5*sin（2*pi*15*t）+2*sin（2*pi*40*t）;%产生时域信号的样本值，向量

Y=fft（x,N）;%N点的DFT计算

mag=abs（Y）;%FFT的振幅

phase=unwrap（angle（Y））;%FFT的相位

%1.以下绘制物理频谱图（即正频部分）

fn=（0:

N/2）*fs/N;%频率轴上的离散频率点，起始于0频（对应直流成分），终

%于Nyquist频率fs/2，共N/2+1个频率点

subplot（2,2,1）%将图形窗口分割为2×2的子窗口，并指定第1个子窗口为绘图区

plot（fn,mag（1:

N/2+1））%取出前N/2+1个振幅作图，即正频率分量

xlabel（'频率/Hz'）;ylabel（'振幅'）;

title（'图1：

物理（正频）幅频图'）;gridon%加网格线

%2.以下绘制全频率的幅频图

fn1=（0:

N-1）*fs/N;

subplot（2,2,2）%指定第2个子窗口为绘图区

plot（fn1,mag）;

xlabel（'频率/Hz'）;ylabel（'振幅'）;

title（'图2：

全频率的幅频图'）;gridon

%3.以下绘制正频部分的相频图

subplot（2,2,3）%指定第3个子窗口为绘图区

plot（fn,phase（1:

N/2+1））;

xlabel（'频率/Hz'）;ylabel（'相位'）;title（'图3：

相频图'）;grid

%4.以下移动零频点

Y1=fftshift（Y）;%fftshift移动频率零点，并将Y的左右两部分交换

mag1=abs（Y1）;%重新计算振幅

fn2=fn1-fs/2;%零点移动到fs/2处，故需重新标记频率轴

subplot（2,2,4）;%指定第4个子窗口为绘图区，最终4幅图绘制在一张图上了

plot（fn2,mag1）;

xlabel（'频率/Hz'）;ylabel（'振幅'）;

title（'图4：

fftshift后的幅频图'）;grid

运行结果如下：

图说明：

1是物理谱图（正频部分），从中看到，该信号包含两个频率15Hz和40Hz。

由于使用的采样频率fs=100Hz，所以Nyquist频率为50Hz，在图2中明显能看到整个频谱图关于Nyquist频率对称，不过Nyquist频率右边的谱图实际上是负频部分，没有意义。

图4是fftshift之后的幅频图，由于它是图2结果的左右交换，因此图2右边变成了负频。

另外，图中的振幅不是真实的信号振幅，从信号x（t）的表达式我们知道15Hz和40Hz这两种频率成分的振幅分别是0.5和2。

要得到真实的振幅，只需要将程序中的mag除以N/2即可。

（3）ifft函数

执行离散傅里叶变换的逆变换，格式

x=ifft（Y）或者x=ifft（Y,m）

Y是FFT的输出结果，返回值x是时域上的结果，m仍然是计算使用的数据点数。

在上例中若程序末尾使用:

xx=ifft（Y,N），则得到采样时刻点上，信号x（t）的样本值。

三、实验内容及要求

实验项目：

给定采样频率51.2Hz及采样点数N=512，计算矩形函数

的振幅频谱，并与理论计算结果对比。

A.显然该信号x（t）是无限长的非周期信号，因此做FFT计算时必须先将信号x（t）截断为有限长度。

令采样频率为fs，采样点数目N，则截断长度是：

T0=N/fs（即总的采样时间）

因此截断长度和采样点数目N成正比。

对于无限长的非周期信号，截断长度应尽可能的大，以接近实际信号，避免结果失真；如果是周期信号，则要求截断长度为信号周期的整数倍，以免出现频谱的“泄漏”。

若给定采样点数N=512，则时间采样序列可用向量t表示：

t=（0:

N-1）/fs，矩形函数x（t）的样本值可以使用MATLAB提供的符号函数sign（请使用helpsign命令查询sign函数的定义）来表示：

x=0.5-0.5*sign（t-1），然后使用fft命令即可获得DFT计算结果。

B.根据

（1）式不难算出上述信号x（t）的傅里叶变换的理论结果（精确值）：

其振幅为：

然后再根据（8）式即可绘出x（t）的振幅频谱，将其与fft计算结果比较：

改变N以改变截断长度观察FFT结果与（8）式结果的差异。

要求：

⑴认真阅读例1提供的程序；

⑵参考例1中的程序，编写出本实验项目的fft计算程序，绘制出零点移动到Nyquist频率处的振幅频谱1；对图形进行标注（如参考图所示）2；使用（8）式绘制矩形函数x（t）精确的振幅频谱图3；使用subplot命令将2幅图在一张图上显示4；频率区间[-5,5]观察比较得到的两个振幅频谱5；改变采样数N，再比较两幅图的差异，分析采样点数对fft计算结果精确度的影响6。

附实验结果参考图：

边瓣

主瓣