实验二 离散时间傅里叶变换

数字信号处理实验

实验二离散时间傅里叶变换

一、实验目的

理解数值计算在离散时间傅里叶变换（DTFT）中的作用。

二、实验原理

经由正、逆离散时间傅里叶变换表达的信号傅里叶表示式是信号分析的一个关键部分。下面方程分别是分析方程与综合方程：

X(e)

j

n

x[n]e

j n

(3.9)

x[n]

12

X(ej )ej nd

(3.10)

类似的，当LTI系统用于滤波的时候，作为冲激响应离散时间傅里叶变换的频率响应，提供了LTI系统简洁的描述。离散时间傅里叶变换X(e)是 的周期复值函数，周期总是

j

2 ，并且基周期通常选在区间[ , )上。

应用MATLAB时需注意： 1、DTFT的定义对无限长信号有效，但当能从变换定义式推导出解析式并只是计算它时，则可以使用MATLAB计算；

2、MATLAB擅长在有线网格点上计算DTFT，通常在[ , )上选择一组均匀隔开的频率，或者对共轭对称变换选择[0, ]区间，这样(3.9)式变为

j k

j2 k/N

X(e) X(e

) x[n]e j(2 k/N)n,k 0,1, N 1

n 0

L 1

(3.11)

DTFT的周期性意味着在 0区间上的数值是对那些k N/2的数值。因为上式是在有限数量的频率点 k 2 k/N处计算，并在有限区间内求和，因此它是可计算的。由于信号长度必须是有限的(0 n L)，这个求和式不适用于x[n] au[n]的情形。若不考虑DFT的物理性质，而只着重于计算DTFT的样本，可计算的上式恰是N点DTFT的离散傅里叶变换（DFT）。 在正确应用FFT计算N点DFT之前，需要对x[n]进行时间混叠，因而此时应确保N L。

n

三、实验内容

1、脉冲信号的DTFT

设矩形脉冲r[n]由下式定义

1,r[n]

0

0 n L其他

（3.12）

a.证明r[n]的DTFT可由下面的数学表示式得出

数字信号处理实验

1sin( L)

R(ej ) e j (L 1)/2sin( )

21sin( L)

asinc( ,L)

sin( )

2

(3.13)

该变换的第一项时常与DTFT具有相关的特殊形式，成为混叠sinc函数：

(3.14)

b.使用dtft函数计算12点脉冲信号的DTFT，绘出在区间 上对 的DTFT，把实部和虚部分开绘出。另绘出DTFT的幅值。选择频率样本的数量是脉冲长度的5到10倍，以使绘出的图看上去很平滑。用不同的频率样本做实验。

c.注意asinc函数的零点位置是规则分布的。对奇数长脉冲，如L=15的脉冲重复进行DTFT计算并绘出幅度同样再次检验零点位置，注意峰值高度。

d.对asinc函数零点的间距与asinc函数的直流值，确定出通用规则。 a.证明如下：

DTFT(r[n])

j L/2

j L/2

n

r[n]e

j L/2

j n

r[n]e

n 0

L 1

j n

1 e j L

1 e j

ee e

e j /2ej /2 e j /2

1 sin( L) e j (L 1)/2 R(ej )sin( )

2

因此r[n]的DTFT可以由式（3.13）得出。

b.DTFT的计算可用下面的函数： M文件：

function [H,W]=dtft(h,N)

%DTFT calculate DTFT at N equally spaced frequencies %Usage:

%[H,W]=dtft(h,N)

%h: finite-length input vector,whose length is L

%N: number of frequencies for evaluation over [pi,pi) %=>constraint N>=L

%H:DTFT values(complex)

%W:(2nd output)vector of freqs where DTFT is computed N=fix(N);

L=length(h); h=h(:); % for vector ONLY if(N<L)

error('DTFT:#Ddata samples cannot exceed # freq samples') end

W=(2*pi/N)*[0:(N-1)]'; mid=ceil(N/2)+1;

W(mid:N)=W(mid:N)-2*pi; %move[pi,2pi) to [-pi,0)

数字信号处理实验

W=fftshift(W);

H=fftshift(fft(h,N)); %move negative freq components

要计算12点脉冲信号的DTFT,则取L=12，并依此取不同的N值，N=60,72,84,96,108,120， 当N=60时，在命令窗口输入: >> r=ones(12,1); >> [H,W]=dtft(r,60); >>subplot(311) >>plot(W,real(H));

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313) >>plot(W,abs(H));

>>grid,xlabel('w'),ylabel('幅度') 得到结果为：

1510

real(H)

50-5-4105

-3

-2

-1

0w

1

2

3

4

imag(H)

0-5-10-41510

-3

-2

-1

0w

1

2

3

4

幅度

50-4

-3-2-1

0w

1234

图2-1-b N=60

当N=72时，

>> [H,W]=dtft(r,72); >>subplot(311)

数字信号处理实验

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313) >>plot(W,abs(H));

>>grid,xlabel('w'),ylabel('幅度') 得到结果为：

real(H)

-4-3

-2

-1

0w

1

2

3

4

imag(H)

-4-3

-2

-1

0w

1

2

3

4

幅度

-4

-3-2-1

0w

123

4

图2-1-b N=72

当N=84时，

>> [H,W]=dtft(r,84); >>subplot(311) >>plot(W,real(H));

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313)

数字信号处理实验

>>grid,xlabel('w'),ylabel('幅度') 得到的结果为：

real(H)

-4-3

-2

-1

0w

1

2

3

4

imag(H)

-4-3

-2

-1

0w

1

2

3

4

幅度

-4

-3-2-1

0w

123

4

图2-1-b N=84

当N=96时，

>> [H,W]=dtft(r,96); >>subplot(311) >>plot(W,real(H));

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313) >>plot(W,abs(H));

>>grid,xlabel('w'),ylabel('幅度') 得到的结果为：

数字信号处理实验

1510

real(H)

50-5-4105

-3

-2

-1

0w

1

2

3

4

imag(H)

0-5-10-41510

-3

-2

-1

0w

1

2

3

4

幅度

50-4

-3-2-1

0w

1234

图2-1-b N=96

当N=108时，

>> [H,W]=dtft(r,108); >>subplot(311) >>plot(W,real(H));

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313) >>plot(W,abs(H));

>>grid,xlabel('w'),ylabel('幅度') 得到的结果为：

数字信号处理实验

real(H)

-4-3

-2

-1

0w

1

2

3

4

imag(H)

-4-3

-2

-1

0w

1

2

3

4

幅度

-4

-3-2-1

0w

123

4

图2-1-b N=108

当N=120时，

>> [H,W]=dtft(r,120); >>subplot(311) >>plot(W,real(H));

>>grid,xlabel('w'),ylabel('real(H)') >>subplot(312) >>plot(W,imag(H));

>>grid,xlabel('w'),ylabel('imag(H)') >>subplot(313) >>plot(W,abs(H));

>>grid,xlabel('w'),ylabel('幅度')

得到的结果为：

数字信号处理实验

1510

real(H)

50-5-4105

-3

-2

-1

0w

1

2

3

4

imag(H)

0-5-10-41510

-3

-2

-1

0w

1

2

3

4

幅度

50-4

-3-2-1

0w

1234

图2-1-b N=120

结果分析：由以上几个不同N值下得到的幅频响应可知，N越大，幅频响应曲线越平滑。 c.分别取不同的L值进行检验，L=15，17,21 当L=15时， >> r=ones(15,1); >> [H,W]=dtft(r,120); >>plot(W/2/pi,abs(H));

>>grid,xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|') 得到的结果为：

数字信号处理实验

15

10

|H(w)|

50-0.5

-0.4-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

0.30.40.5

图2-1-c L=15

当L=17时， >> r=ones(17,1); >> [H,W]=dtft(r,200); >>plot(W/2/pi,abs(H));

>>grid,xlabel('NORMALIZED FREQUENC'),ylabel('幅度') 得到的结果为：

1816

141210

幅度

86420-0.5

-0.4-0.3

-0.2-0.100.10.2

NORMALIZED FREQUENC

0.30.40.5

图2-1-c L=17

当L=21时， >> r=ones(21,1); >> [H,W]=dtft(r,200);

数字信号处理实验

>>plot(W/2/pi,abs(H));

>>grid,xlabel('NORMALIZED FREQUENC'),ylabel('幅度')

得到的结果为：

25

20

15

幅度

1050-0.5

-0.4-0.3

-0.2-0.100.10.2

NORMALIZED FREQUENC

0.30.40.5

图2-1-c L=21

结果分析：（1）由上三图可知，当L为奇数时，零点间距是固定的。 （2）而其峰值即为L的大小。

d.根据以上分析，asinc函数零点的间距为 2 /L，峰值为L。

2.、asinc函数的M文件

编写一个MATLAB的函数如asinc（W，L），直接从（3.6）式计算在频率格上的asinc（W，L）。该函数应有两个输入：长度L和频率W的向量。函数必须检查被零除的情形，如w=0时。

直接计算混叠sinc函数（3.14）式得到脉冲信号的DTFT。绘出幅度，保存该图以便与dtft得到的结果比较。 编写的asinc函数如下： function a=asinc(w,L) if (w==0) a=L;

else a=sin(1/2*w*L)/sin(1/2*w); end

利用混叠sinc函数式（3.14）得到DTFT，M文件： function R=a(w,L)

% w:数字域频率基本频率在[-pi,pi]之间 % L：矩形脉冲的长度

b=sin(1/2*w*L)./sin(1/2*w); % 混叠sinc函数 c=exp(-j*w*(L-1)/2);

数字信号处理实验

R=b.*c; % 矩形脉冲的DTFT subplot(211)

plot(w/2/pi,abs(R)) % 频谱 grid

xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|'),title('由数学式得出的幅频特性') subplot(212)

plot(w/2/pi,180/pi*angle(R)) %相频特性 grid

xlabel('NORMALIZED FREQUENCY'),ylabel('DEGREES'),title('由数学式得出的相频特性') 取L=5，在命令窗口输入： >> w=[-pi:pi/50:pi]; >> L=5; >>a(w,L); 得到的结果为：

由数学式得出的幅频特性

64

20-0.5

|H(w)|

-0.4-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY由数学式得出的相频特性

0.30.40.5

200

DEGREES

100

0-100-200-0.5

-0.4

-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

0.3

0.4

0.5

图2-2 由数学式得到的结果

通过dtft得到DTFT： >> n=0:4;

>> r=ones(5,1);

>> [H,W]=dtft(r,128); >>subplot(211)

>>plot(W/2/pi,abs(H));

>>grid,title('幅频响应'),xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|') >>subplot(212)

>>plot(W/2/pi,180/pi*angle(H));

>>grid,title('相频响应'),xlabel('NORMALIZED FREQUENC'),ylabel('DEGREES')

得到的结果为：

数字信号处理实验

幅频响应

-0.5

|H(w)|

-0.4-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

相频响应

0.30.40.5

DEGREES

-0.5

-0.4

-0.3

-0.2-0.100.10.2NORMALIZED FREQUENC

0.3

0.4

0.5

图2-2 由dtft得到的结果

结果分析：通过比较上面两个图，可知用两种方法得到的DTFT的结果相同。

3、无限长信号的DTFT

通常不可能计算一个无限长信号的DTFT，但指数信号h[n] au[n]计算比较容易。 当|a|<1, 有

n

h[n] anu[n]

H(e) anu[n]e j n

j

n 0

11 ae j

(3.7)

利用freqz函数可以计算上式：

[HH,WW] freqz(b,a,N,'whole')

与dtft类似，freqz有两个输出：交换数值（HH）和频率格点（WW），第四个输入参数

是可以选择的，但如将其设定为’whole’,则输出向量WW指定频率格点的范围是从 0到2 。如果省略第四个参数，频率格点由 区间上等间距的N个点组成。

4、指数信号

对于信号x[n] 0.9u[n]，使用freqz函数计算其DTFTX(e)。

a. 对 在区间 上绘出幅度与相位特性。这需要从freqz函数返回的[X,W]

向量的移位。解释为什么幅度特性是 的偶函数，而相位特性是 的奇函数。

b.推算一阶系统的幅度特性与相位特性的表达式。

c.直接以这些表达式来计算幅度特性与相位特性，并与用freqz函数计算出的结果相

对比。 a.程序：

n

j

数字信号处理实验

>> a=[1,-0.9]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> WW=fftshift(WW); >> HH=fftshift(HH);

>>subplot(211),plot(WW/2/pi,abs(HH));

>>grid,title('幅频响应'),xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|') >>subplot(212),plot(WW/2/pi,180/pi*angle(HH));

>>grid,title('相频响应'),xlabel('NORMALIZED FREQUENCY'),ylabel('DEGREES') 结果：

幅频响应

1510

50-0.5

|H(w)|

-0.4-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

相频响应

0.30.40.5

100

DEGREES

50

0-50-100-0.5

-0.4

-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

0.3

0.4

0.5

图2-4-a

结果分析：由上图可知，幅度特性是 的偶函数，相位特性是 的奇函数。

因为序列x(n)可以分解为共轭对称分量xe(n)和共轭反对称分量xo(n)之和，即

j j

x(n) xe(n) xo(n)，则X(ej ） X（ jX（，对于实序列，只有Re）Ie）

j X（e）R

，因此

j *-j

X(ej ） X（e） X（e）R

，即

j -j j -j

X（ X（），X（ X（），因此幅度是 的偶函数，相位Re）ReIe）Ie

是 的奇函数。

b.该一阶系统的系统函数为：H(ej )

1j jarg[H(ej )]

，其幅度为 |H(e)|e

1 0.9e j

|H(ei )|，相位为arg[H(ei )]。

数字信号处理实验

c.直接由（b）推出的表达式计算幅度特性和相位特性，程序如下： >> w=-pi:pi/50:pi;

>> y=1./(1-0.9.*exp(-j*w));

>>subplot(211),plot(w/2/pi,abs(y));

>>grid,title('幅频响应'),xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|') >>subplot(212),plot(w/2/pi,180/pi*angle(y));

>>grid,title('相频响应'),xlabel('NORMALIZED FREQUENCY'),ylabel('DEGREES') 结果：

幅频响应

-0.5

|H(w)|

-0.4-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

相频响应

0.30.40.5

DEGREES

-0.5

-0.4

-0.3

-0.2-0.100.10.2NORMALIZED FREQUENCY

0.3

0.4

0.5

图2-4-c

结果分析：由b、c得到的图形可知，结果相同。

5、复指数信号

如果取（3.7）中a z0 rej 为复数，则可应用相同的变换。由此可了解这个复数的幅度和相位是如何影响DTFT的。

n

u[n]。用subplot指令将实部和虚部 a.取z0 0.95ej3 /11，对0 n 30，绘出x[n] z0

对n绘出双子图。

j3 /11

b.再一次取z0 0.95e，计算其DTFT，并对 绘出幅度。注意作为 的函数，幅度响

应的尖峰位于何处，把尖峰位置与z0的极坐标形式联系起来。

c.如果z0的角度变为 3 /5，粗略绘出预期的DTFT幅度，并通过对freqz的计算进行绘图来验证。

d.改变幅度r |z0|，重新绘出其DTFT。r取四个数值：r=0.975、0.95、0.9和0.8.注意，

数字信号处理实验

当|z0|接近于1时，DTFT幅度尖峰的高度和带宽都将改变。在-3dB点测量带宽的数值。试建立一个联系带宽与r的简单表达式。 a.程序：

z0=0.95*exp(j*3*pi/11); r=ones(1,31); nn=0:30;

x=(z0.^nn).*r; subplot(211) stem(nn,real(x));

grid,xlabel('n'),ylabel('real(x)') subplot(212) stem(nn,imag(x));

grid,xlabel('n'),ylabel('imag(x)') 结果：

real(x)

0510

15n

202530

imag(x)

0510

15n

202530

图2-5-a

b.程序：

>> a=[1,-0.9*exp(j*3*pi/11)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> WW=fftshift(WW); >> HH=fftshift(HH);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

数字信号处理实验

幅频响应

|H(w)|

-4

-3

-2

-1

0w

1

2

3

4

图2-5-b

结果分析：由上图可知尖峰出现的位置为0.85左右，而3 /11 0.8567。可见尖峰出现的位置即为其极坐标的对应的角度。 c.改变角度后， 程序：

>> a=[1,-0.9*exp(j*3*pi/5)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> WW=fftshift(WW); >> HH=fftshift(HH);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

数字信号处理实验

幅频响应

|H(w)|

-4

-3-2-1

0w

1234

图2-5-c

结果分析：由上图可知，在w=1.88处出现尖峰，与预期的DTFT幅度大致一样。 d.改变幅度， 当r=0.975时，

>> a=[1,-0.975*exp(j*3*pi/5)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> WW=fftshift(WW); >> HH=fftshift(HH);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

数字信号处理实验

幅频响应

-4

|H(w)|

-3-2-1

0w

1234

图2-5-d r=0.975

当r=0.95时，

>> a=[1,-0.95*exp(j*3*pi/5)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> HH=fftshift(HH); >> WW=fftshift(WW);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

数字信号处理实验

幅频响应

2018

161412

|H(w)|

1086420-4

-3

-2

-1

0w

1

2

3

4

图2-5-d r=0.95

当r=0.9时，

>> a=[1,-0.9*exp(j*3*pi/5)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> HH=fftshift(HH); >> WW=fftshift(WW);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

幅频响应

12

10

8

|H(w)|

6

4

2

0-4

-3-2-1

0w

1234

图2-5-d r=0.9

数字信号处理实验

当r=0.8时，

>> a=[1,-0.8*exp(j*3*pi/5)]; >> b=1; >> N=100;

>> [HH,WW]=freqz(b,a,N,'whole'); >>mid=ceil(N/2)+1;

>>WW(mid:N)=WW(mid:N)-2*pi; >> HH=fftshift(HH); >> WW=fftshift(WW);

>>subplot(111),plot(WW,abs(HH));

>>grid,title('幅频响应'),xlabel('w'),ylabel('|H(w)|') 结果：

幅频响应

5.55

4.543.5

|H(w)|

32.521.510.5-4

-3

-2

-1

0w

1

2

3

4

图2-5-d r=0.8

结果分析：由上图可知，当|z0|接近于1时，DTFT幅度尖峰的高度和带宽都将改变。

本文来源：https://www.bwwdw.com/article/2ot1.html