实验用MATLAB计算傅里叶变换

实验二 用MATLAB 计算傅立叶变换
（2课时）
一、实验目的
1、掌握用MA TLAB 计算DTFT 及系统频率响应的方法。

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

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

二、实验设备
计算机一台，装有MATLAB 软件。

三、实验原理和基本操作
1．用MA TLAB 计算DTFT
对于序列x （n ），其离散时间傅立叶变换（DTFT ）定义为：
∑∞-∞=-=
n n j e n x j X ωω)()( (1)
序列的傅立叶变换（DTFT ）在频域是连续的，并且以ω=2π为周期。

因此只需要知道jw X(e )的一个周期，即ω=[0，2π]，或[-π，π]。

就可以分析序列的频谱。

用MA TLAB 计算DTFT ，必须在-π≤ω≤π范围内，把ω用很密的、长度很长的向量来近似，该向量中各个值可用下式表示： w=k*dw=k*
K π2 （2） 其中：d ω=K
π2 称为频率分辨率。

它表示把数字频率的范围2π均分成K 份后，每一份的大小，k 是表示频率序数的整数向量，简称为频序向量，它的取值可以有几种方法：通常在DTFT 中，频率取-π≤ω<л的范围，当K 为偶数时，取 k 12,,1,0,1,,12,2--+--
=K K K 如果K 为奇数，则取 k 5.02
,,1,0,1,,5.02--+-=K K 可以为奇偶两种情况综合出一个共同的确定频序向量k 的公式； k=12K -⎢⎥-⎢⎥⎣⎦ ：12K -⎢⎥⎢⎥⎣⎦
（3） 上式中⎢⎥⎣⎦表示向下取整。

在MA TLAB 中的向下取整函数为floor ，floor （x ）的作用是把x 向下（向-∞方向）取整，所以与（3）式等价的MATLAB 语句为 k ))5.02
(:)5.02((-+-=K K floor （4） 给定了输入序列（包括序列x 及其位置向量n ），又设定了频率分辨率d ω及频序向量k ，则DTFT 的计算式（1）可以用一个向量与矩阵相乘的运算来实现。

⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡=---------N K N N K K n j n j n j n j n j n j n j n j n j N K e e e e e e e e e n x n x X X X ωωωωωωωωωωωω 212222111211)](,),([)](,),(),([121 （5） 如果频率向量表为ω=[ω1,ω2,……，ωκ]=k*d ω,而序列的位置向量为nx=[n1：nN]，则
（5）式中的矩阵的指数部分可以写成-j*n T *ω，用MA TLAB 语句表达时，把ω代以w ，转置符号n T 换成MA TLAB 中的相应符号n '
，则求DTFT 的程序可以写成： )'***exp(*k n dw j x X -=
例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． 用MA TLA B 计算系统频率响应()jw H e
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≤ω<π之间的特性，如果要让MATLAB 计算全频率0≤ω<2π范围的频率特性，要增加一个输入变元’whole ’： [H ，w]=freqz(b,a,N,’whole ’),此时w 的范围是0到2π（不含2π）。

若没有左端变量，即键入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 )=ωω
j j e
e ---+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． 用MA TLAB 计算DFT 和IDFT
可以利用MATLAB 提供的计算快速傅立叶变换（FFT ）和快速傅立叶反变换（IFFT ）的函数来计算离散傅立叶变换（DFT ）和离散傅立叶反变换（IDFT ）。

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

X=fft(x)完成对序列x 的L 点DFT ，其中L 为序列x 的长度。

X=fft(x ，N)则指定了采用N 点DFT 。

如果N>L,则程序会自动给x 后面补零,使其长度为N;如果N<L,则程序会自动将x 截断,取前N 个数据。

同样，计算离散傅立叶反变换时，调用函数ifft 方法为ifft （X ）或ifft(X ，N)，其输入变元的意义与fft 函数相同。

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

解：x=[2,-1,1,1];
X=fft(x)
运行结果：
4．用MA TLAB 计算圆周卷积
圆周卷积定理：设有限长序列x 1(n)和x 2(n)的长度分别为N 1和N 2，N=max[N 1, N 2]。

x 1(n)和x 2(n)的N 点DFT 分别为
X 1（k ）=DFT[x 1(n)], X 2(k)=DFT[x 2(n)].
若X （k ）= X 1（k ）* X 2(k)
则x 1(n)与x 2(n)的N 点圆周卷积是
x(n)=IDFT[X(k)]=
)())(()(2101n R m n x m x N N N m -∑-=
可以利用此定理和函数fft 、ifft 计算圆周卷积
例4：设x 1(n)={1，2，3}，x 2(n)={5，4，-3，-2}，计算4点圆周卷积y(n)= x 1(n)④ x 2(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．利用圆周卷积计算线性卷积
设x 1(n)为N 1点序列，x 2(n)为N 2点序列，x 1(n)和x 2(n)的线性卷积为
y(n)=
)()(21
011m n x m x N m -∑-=
y(n)为（N 1+ N 2-1）点序列。

如果N ≥ N 1+ N 2-1，则N 点圆周卷积能代表线性卷积。

例5：设x 1(n)与x 2(n)是两个四点序列：
x 1(n)={1，2，2，1}， x 2(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程序为
n=0:9;x=n+1;h=[1,0,-1];N=6;
y1=ovrlpsav(x,h,N)
y2=fftfilt(h,x)
运行结果：
四、实验内容
1、上机运行上面给出的各个例子，记录程序运行结果。

2、求有限长序列x(n)=[1,3,5,2,1]的DTFT，画出它在ω= -8~8rad/s范围内的幅频和相频特性曲线。

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

五、注意事项
1、实验前复习相关理论知识。

2、注意子函数调用时输入输出参数的匹配。