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

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

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 )ttt011(3)s f f T NT N∆===因此要增加分辨率须增加采样点数目N 。

频域上每个离散点对应的频率为:0; 0,1,2,...,-1 (4)s n f n nf n n N T NT N====显然n = 0的点对应于直流成分。

经过以上离散化处理之后，连续积分的傅里叶变换(1)式转变为如下离散形式：12/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)式通常又简写成如下形式：10()(), 0,1,2,..., 1 (6)N nkN 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 的思想实质是fff()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)Ｙ返回值（复数），返回m点的DFT序列，即(6)式左边的X；m计算时使用的数据点数（样本数）；x时域信号x(t)在采样点t k处的值，即(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 的左右两半。

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

例１：对时域信号()0.5sin(215)2sin(240)x t t t ππ=⋅+⋅进行频谱分析。

fs=100; % 采样频率>2倍的信号频率 N=256; % 采样点数目（=2的幂次） n=0:N-1; % 构造采样点序列t=n/fs; % 得到采样时间序列，t=nT=n/fsx=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个振幅作图，即正频率分量……f n ： f 轴1T 02T 0/22s N f T =1N T -s Nf T = 对称轴(Nyquist 频率)负频部分→← 正频部分 FFT ：注意：频率轴上最后一个点的频率等于采样频率， 实际上不存在。

0, (0,1,2,...,1)n s n n n f f n N T NT N====- 为频率轴上的频率点。

特别： 2212N s s N f f f N == 为Nyquist 频率，整个FFT 频谱关于此频率对称， 右边实际上是负频xlabel('频率/Hz');ylabel('振幅');title('图1：物理(正频)幅频图');grid on % 加网格线% 2. 以下绘制全频率的幅频图fn1=(0:N-1)*fs/N;subplot(2,2,2) %指定第2个子窗口为绘图区plot(fn1,mag);xlabel('频率/Hz');ylabel('振幅');title('图2：全频率的幅频图'); grid on%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 。