实验用MATLAB计算傅里叶变换

实验二 用MATLAB计算傅立叶变换

（2课时）

一、实验目的

1、掌握用MATLAB计算DTFT及系统频率响应的方法。

2、掌握用MATLAB计算DFT和IDFT的方法。

3、掌握用DFT计算圆周卷积和线性卷积的方法。

二、实验设备

计算机一台，装有MATLAB软件。

三、实验原理和基本操作

1．用MATLAB计算DTFT

对于序列x（n），其离散时间傅立叶变换（DTFT）定义为：

nnjenxjX)()( (1)

序列的傅立叶变换（DTFT）在频域是连续的，并且以=2π为周期。因此只需要知道jwX(e)的一个周期，即=[0，2π]，或[-π，π]。就可以分析序列的频谱。

用MATLAB计算DTFT，必须在-π≤≤π范围内，把用很密的、长度很长的向量来近似，该向量中各个值可用下式表示：

w=k*dw=k*K2 （2）

其中：d=K2 称为频率分辨率。它表示把数字频率的范围2π均分成K份后，每一份的大小，k是表示频率序数的整数向量，简称为频序向量，它的取值可以有几种方法：通常在DTFT中，频率取-π≤

k12,,1,0,1,,12,2KKK

如果K为奇数，则取

k5.02,,1,0,1,,5.02KK

可以为奇偶两种情况综合出一个共同的确定频序向量k的公式；

k=12K ：12K （3）

上式中表示向下取整。 在MATLAB中的向下取整函数为floor，floor（x）的作用是把x向下（向-方向）取整，所以与（3）式等价的MATLAB语句为

k))5.02(:)5.02((KKfloor （4）

给定了输入序列（包括序列x及其位置向量n），又设定了频率分辨率d及频序向量k，则DTFT的计算式（1）可以用一个向量与矩阵相乘的运算来实现。

NKNNKKnjnjnjnjnjnjnjnjnjNKeeeeeeeeenxnxXXX212222111211)](,),([)](,),(),([121 （5）

如果频率向量表为=[1,2,……，]=k*d,而序列的位置向量为nx=[n1：nN]，则（5）式中的矩阵的指数部分可以写成-j*nT*，用MATLAB语句表达时，把代以w，转置符号nT换成MATLAB中的相应符号n'，则求DTFT的程序可以写成：

)'***exp(*kndwjxX

例1 求有限序列x=[2，-1，1，1]的DTFT，其位置向量为nx=[0：3]。

假如取64个频点，画出它在-ππ范围内的幅频和相频特性。程序如下：

x=[2,-1,1,1];nx=[0:3];K=64;dw=2*pi/K;

k=floor((-K/2+0.5):(K/2-0.5)); % 设定频序向量

%w=linspace(-8,8,1000);

X=x*exp(-j*dw*nx'*k) % 用（1）式计算DTFT

subplot(2,1,1); plot(k*dw,abs(X)), % 画幅频曲线

subplot(2,1,2) ;plot(k*dw,angle(X)), % 画相频曲线

2． 用MATLA B计算系统频率响应()jwHe

MATLA B中用于求系统频率响应的函数是freqz,其功能是：由给定的系统函数H(Z)的分子和分母的系数向量绘制系统的幅度和相位响应。

调用格式：

[H，w]=freqz(b,a,N)，或N缺省[H，w]=freqz(b,a),此时N取默认值512。

H是系统的频率特性，它是一个N元的复数向量。w是数字频率向量，它把0到π均分为N份，分辨率为π/N，w =[0：N-1]π/N。

b和a分别为分子分母多项式的负幂系数向量，即多项式的首项是常数项，以后按ej

的升幂排列，由此形成的多项式的系数向量。 N为所选的频率点数，它决定了频率分辨率的密度。

这样求出的频率特性是在0≤

若没有左端变量，即键入freqz(b,a,N)，MATLAB将不给出数据H和w而只绘出频率特性曲线。若在输入变元中给出频点向量w: freqz(b,a,w)，这种调用方法可以在自己选定的频点向量w上计算频率特性。

例2：设一个LTI系统的差分方程为y(n)-0.9y(n-1)=0.5x(n)+0.8x(n-1)

求其频率响应。

可直接写出H(ej)=jjee9.018.05.0

用下面的程序可绘出幅频和相频特性曲线。

b=[0.5 0.8]; % 分子多项式系数向量

a=[1 -0.9]; % 分母多项式系数向量

[H,w]=freqz(b,a); % 求出频率响应（0到pi分成500点）

subplot(2,1,1),plot(w,abs(H)),grid on

xlabel('频率'),ylabel('幅度dB')

subplot(2,1,2),plot(w,angle(H)),grid on

xlabel('频率'),ylabel('相角(度)')

3． 用MATLAB计算DFT和IDFT

可以利用MATLAB提供的计算快速傅立叶变换（FFT）和快速傅立叶反变换（IFFT）的函数来计算离散傅立叶变换（DFT）和离散傅立叶反变换（IDFT）。

在MATLAB信号处理工具箱中，函数fft和ifft用于快速傅立叶变换和反变换。

X=fft(x)完成对序列x的L点DFT，其中L为序列x的长度。X=fft(x，N)则指定了采用N点DFT。如果N>L,则程序会自动给x后面补零,使其长度为N;如果N

例3：长度为4的有限序列：x(0)=2,x(1)= -1,x(2)=1,x(3)=1,求它的DFT。

解：x=[2,-1,1,1];

X=fft(x)

运行结果： 4．用MATLAB计算圆周卷积

圆周卷积定理：设有限长序列x1(n)和x2(n)的长度分别为N1和N2，N=max[N1, N2]。x1(n)和x2(n)的N点DFT分别为

X1（k）=DFT[x1(n)], X2(k)=DFT[x2(n)].

若X（k）= X1（k）* X2(k)

则x1(n)与x2(n)的N点圆周卷积是

x(n)=IDFT[X(k)]=)())(()(2101nRmnxmxNNNm

可以利用此定理和函数fft、ifft计算圆周卷积

例4：设x1(n)={1，2，3}，x2(n)={5，4，-3，-2}，计算4点圆周卷积y(n)= x1(n)④ x2(n)

解：写出MATLAB程序如下：

x1=[1 2 3],x2=[5 4 -3 -2]

X1 =fft(x1 ,4); X2 =fft(x2 ,4);Y= X1 .* X2 ;

y=ifft(Y,4)

运行结果：

5．利用圆周卷积计算线性卷积

设x1(n)为N1点序列，x2(n)为N2点序列，x1(n)和x2(n)的线性卷积为

y(n)= )()(21011mnxmxNm

y(n)为（N1+ N2-1）点序列。

如果N N1+ N2-1，则N点圆周卷积能代表线性卷积。

例5：设x1(n)与x2(n)是两个四点序列：

x1(n)={1，2，2，1}， x2(n)={1，-1，1，-1}

(1) 求它们的线性卷积y(n);

(2) 计算圆周卷积，使得它与y(n)相等。

解：线性卷积可以调用conv函数来求。圆周卷积可按例4的方法求。程序如下：

x1=[1,2,2,1];x2 =[1,-1,1,-1];

y=conv(x1 , x2 )

X1 =fft(x1 ,7); X2 =fft(x2 ,7); Y= X1 .* X2 ; x3 =ifft(Y,7)

运行结果：

6．分段卷积

当x(n)长度远远大于h(n)时，要采用分段卷积的方法，以减少线性卷积的运算量。分段卷积的方法有重叠保留法和重叠相加法两种。

(1) 重叠保留法

设h(n)的长度为M，x(n)的长度为Lx，分段后每一段的有效数据的长度为L=N-M+1。编写出实现重叠保留法的函数ovrlpsav如下：

function [y]=ovrlpsav(x,h,N)

Lx=length(x);M=length(h);

M1=M-1;L=N-M1;

H=fft(h,N);

x=[zeros(1,M1),x,zeros(1,N-1)];

K=floor((Lx+M1-1)/(L))+1;

Y=zeros(K+1,N);

for k=0:K-1

xk=x(k*L+1:k*L+N);

Xk=fft(xk);

Y(k+1,:)=real(ifft(Xk.*H));

end

Y=Y(:,M:N)';

y=(Y(:))';

(2) 重叠相加法

MATLAB信号处理工具箱中的函数fftfilt.m可实现重叠相加法。它有两种调用格式：

1）y=fftfilt(h,x)

2）y=fftfilt(h,x,r)

其中h代表系统的脉冲响应h(n),x是输入序列。在第一种格式中，程序自动把输入分成每段512个样本。并按512点（如果h(n) 的长度比512长，则按h(n)的长度）FFT进行各段的卷积运算。在第二种调用格式中，r是用户指定的FFT长度，而输入x就按这个长度分段。

例6：设x(n)=n+1,0n9,h(n)={1,0,-1}。分别用ovrlpsav(按N=6)和fftfilt函数求x(n)与h(n)的线性卷积

解：MATLAB程序为