哈尔滨工程大学 数字信号处理实验二

实 验 报 告

课程名称 实验项目名称 实验类型 班级 姓名 实验室名称 预习部分 实验成绩 教师签字

实验过程 表现 实验学时 学号 指导教师 实验时间 实验报告 部分 日期 总成绩

实验二 离散时间傅立叶变换 实验报告

一、实验原理

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

X(ej?)=

n????x[n]eπ??j?n

1j?j?nx[n]?X(e)ed? ?2π?π类似地，当LTI系统用于滤波时，作为冲击响应离散时间傅里叶变换的频率响应，提供了LTI系统简介的描述。离散时间傅里叶变换X(ej?)是?的周期复值函数，周期总是

2π，并且基周期通常选在区间[-π，π）上。对离散时间傅里叶变换DTFT来说有两个问题：

第一个问题是：DTFT的定义对无限长信号是有效的 第二个问题是：DTFT是连续变量?的函数

在MATLAB中，任何信号必须是有限长度的，这就使第一点成为问题。因此，不可能使用MATLAB计算无限长信号的DTFT。有一个值得注意的例外情形，当能从变换定义式推导出解析式并只是计算它时，可以使用MATLAB计算无限长信号的DTFT。例如在x[n]=au[n]，x[n]具有有理的DTFT的情形下。

2、第二个问题是频率抽样问题。MATLAB擅长在有限网格点上计算DTFT。通常选择足够多的频率以使绘出的图平滑，逼近真实的DTFT。对计算有利的最好选择是在（-π，π）区间上一组均匀地隔开的频率，或者对共轭对称变换选择[0,π]区间。采用上述抽样办法，DTFT式变成X(ej?n)=X(ej2πk/N)??x[n]e?j(2πk/N)n,k?0,1,2...N?1

n?0L?1（3.11）DTFT的周期性意味着在-π≤?<0区间上的数值是那些对k>N/2的数值。因为上式是在有限数量的频率点?k=2πk/N处计算，并在有限范围内求和，因此它是可计算的。由于信号长度必须是有限的（0≤n

算法）时，如果N=L很便利。的确，在使用N点DFT时，如果N>L，只须想像x[n]是增补零的即可；如果N=L）点数多得多的频率处来计算DTFT。

3、计算DTFT需要两个函数，MATLAB的freqz函数计算无限长信号，dtft（h，H）函数计算有限长信号的DTFT。

n

二：实验内容

实验内容包括脉冲信号的DTFT，asinc的M文件，无限长信号的DTFT，指数信号，复指数信号这五个部分。 1. 脉冲信号的DTFT

a、已知矩形脉冲r[n]的定义式为在0≤n

n????x[n]e??j?n中，计算化简后可得：

R(ejw)?{sin(wl/2)/sin(w/2)}*e?jw(L?1)/2。

b、建立DTFT函数（.M文件） 实验程序清单1：

function[A,B]=dtft(h,N)

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

%[A,B]=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 vectors ONLY!!! if(N

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

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

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

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

实验程序清单2:（当w=80时）

h=ones(1,12);

[A,B]=dtft(h,80); %w=80时

subplot(311),plot(B,real(A)); %把绘图窗口分成三行一列三个子图，然后在第一块区域作图 xlabel('B'),ylabel('A,real part'); %以B为横轴，A的实部为纵轴

subplot(312),plot(B,imag(A)); %把绘图窗口分成三行一列三个子图，然后在第二块区域作图 xlabel('B'),ylabel('A,imag part'); %以B为横轴，A的虚部为纵轴

subplot(313),plot(B,abs(A)); %把绘图窗口分成三行一列三个子图，然后在第三块区域作图 xlabel('B'),ylabel('A amplitude'); %以B为横轴，A的幅值为纵轴

实验仿真图：

20A,real part100-10-410-3-2-10B1234A,imag part0-10-415-3-2-10B1234A amplitude1050-4-3-2-10B1234 实验程序清单3:（当w=100时） h=ones(1,12);

[A,B]=dtft(h,100); %w=100时 subplot(311),plot(B,real(A)); xlabel('B'),ylabel('A,real part'); subplot(312),plot(B,imag(A)); xlabel('B'),ylabel('A,imag part'); subplot(313),plot(B,abs(A)); xlabel('B'),ylabel('A amplitude'); 实验仿真图：

20A,real part100-10-410-3-2-10B1234A,imag part0-10-415-3-2-10B1234A amplitude1050-4-3-2-10B1234 实验程序清单4:（当w=120时） h=ones(1,12);

[A,B]=dtft(h,120); %w=120时 subplot(311),plot(B,real(A)); xlabel('B'),ylabel('A,real part'); subplot(312),plot(B,imag(A)); xlabel('B'),ylabel('A,imag part'); subplot(313),plot(B,abs(A)); xlabel('B'),ylabel('A amplitude');

实验仿真图：

20A,real part100-10-410-3-2-10B1234A,imag part0-10-415-3-2-10B1234A amplitude1050-4-3-2-10B1234 实验结果分析：由以上三组不同频率（w=80,w=100,w=120）的仿真图形可知，当w越大，取的点数越多时，绘制的图形越来越平滑。

c、注意asinc函数零点的位置是规则分布的。对奇数长脉冲，比如L=15的脉冲重复进行DTFT计算并绘出幅度；同样再次检验零点位置，注意峰值高度。 实验程序清单： n=0:14;

h=ones(1,15);

[X,W]=dtft(h,150);

subplot(211),plot(W,abs(X)); %把绘图窗口分成二行一列二个子图，然后在第一块区域作图

grid; %图层为网格

xlabel('W'),ylabel('|X(w)|'); %X轴为w，Y轴为X(w)幅值 [H,X]=dtft(X,150);

subplot(212),plot(X,abs(H)); %把绘图窗口分成二行一列二个子图，然后在第二块区域 作图 grid;

xlabel('X'),ylabel('|H(w)|');%X轴为X，Y轴为H(w)幅值

实验仿真图：

15|X(w)|1050-4150-3-2-10W1234100|H(w)|500-4-3-2-1 实验结果分析：由仿真结果可以看出函数零点关于x=0对称。在对称轴处幅值最大，离对称轴越远，幅值逐渐减小。

d、对asinc函数零点的间距与asinc函数的直流值，确定出通用规则。 零点间距为π/6，由数学推导：在a部分中，当L=12时，由R(ej?0X1234)=0可得sin(wL/2)=0，

即wL/2=k* π 则w=k*π/6。故得出以上结论。在实验a部分的仿真图中，可观察到幅值高度为12。则零点间距*幅值=(π/6)*12=2*π。

2、asinc的M文件

编写一个MATLAB的函数如asinc(w,L),直接从混叠函数表达式中计算在频率格上的asinc(w,L)。该函数应该有两个输入：唱的L和频率w的向量。函数必须检查被零除的情况，如w=0时。直接计算混叠sinc函数得到脉冲信号的DTFT。绘出幅度，保存该图以便将其与用dtft得到的结果进行比较。

function [H,W]=asinc(w,L) %定义asinc函数 N=fix(L); if(N<0)

error('Please input a positive integer.(L>0)'); end

%if(sin(w/2)=0) %end

W=(2*pi/N)*[0:(N-1)]'; H=sin(w*N/2)/sin(w/2);

3、指数信号

对于信号x[n]=(0.9)^n*u[n]，使用freqz函数计算其DTFT X(ej?)。

a、对w在区间[-π，π）上绘出幅值与相位特性。这需要从freqz返回的[X,W]向量的移位。

解释为什么幅度特性是w的偶函数，而相位特性是w的奇函数。 实验程序清单1： N=100;

a=[1,-0.9]; b=1;

[H,W]=freqz(b,a,N); W=[-pi:0.1:pi]; H=freqz(b,a,W);

subplot(211);%把绘图窗口分成二行一列二个子图，然后在第一块区域作图 plot(W,abs(H));

grid,title('MAGNITUDE RESPONSE') xlabel('FREQUENCY W'),ylabel('|H(w)|')

subplot(212);%把绘图窗口分成二行一列二个子图，然后在第二块区域作图 plot(W,angle(H));

grid,title('PHASE RESPONSE')

xlabel('FREQUENCY W'),ylabel('DEGREES')

实验仿真图：

MAGNITUDE RESPONSE108|H(w)|6420-4210-1-2-4-3-2-101FREQUENCY WPHASE RESPONSE234DEGREES-3-2-101FREQUENCY W234 j?实验结果分析：将该指数信号代入分析方程X(e)=

n???j??x[n]e??j?n中，计算化简后可从结果

X(e)看出：该结果的幅值关于w是偶函数，相位关于w是奇函数。

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

设一阶差分方程：y[n]=x[n]+ay[n-1],其系统方程H(e其幅度H(e)?j?j?)=

1

1-ae?j?11?a-2acos?2

相位?H(e)??arctan[j?asin?]

1?acos?c、直接以这些表达式来计算幅度特性与相位特性，并与用freqz函数计算出的结果相对比。 x[n]=0.9u[n]的幅度特性为X(ej?)?n1

1.81-1.8cos? 相位特性为?X(e)??arctan[实验程序清单： N=1000; a=0.9;

j?0.9sin?]

1?0.9cos?w=(2*pi/N)*[-(N-1)/2:(N-1)/2]; m=1./((1-a*cos(w))+(a*sin(w).^2)); n=1./((1-a*cos(w)).^2+(a*sin(w)).^2); H=sqrt(m.^2+n.^2);%H的表达式

subplot(2,1,1);%把绘图窗口分成二行一列二个子图，然后在第一块区域作图 plot(w,sqrt(H)); grid;

xlabel(' W'),ylabel('|H(w)|'); phase=(a*sin(w))./(1-a*cos(w)); subplot(2,1,2);

plot(w,atan(-phase));%绘制连续图形 grid,title('PHASE ');

xlabel(' W'),ylabel('DEGREES');

实验仿真图：

1510|H(w)|50-4210-1-2-4-3-2-10 WPHASE 1234DEGREES-3-2-10 W1234 实验结果分析：

以上两个图对比，可以发现用表达式来计算幅度特性和相位特性这种方式与用freqz函数计算这种方式，它们最终得出的结果是一样的。

3、复指数信号：

a、取z=0.95e^j3pi/11,对[0,30]，绘出x[n]=z^n*u[n]。用subplot指令将实部和虚部对

n绘出双子图。 实验程序清单：

z0=0.95*exp(j*3*pi/11); n=[0:30]; x=z0.^n;

subplot(311);%用subplot指令绘出子图

stem(n,x); %以x为横轴，n为纵轴绘制离散信号波形 title('复指数信号');

xlabel('n'),ylabel('x'); subplot(3,1,2); stem(n,real(x)); title('REAL PART');

xlabel('n'),ylabel('real'); subplot(313);

stem(n,imag(x)); %以n为横轴，x的虚部为纵轴绘制离散信号波形 title('IMAG PART');

xlabel('n'),ylabel('imag');

实验仿真图：

复指数信号10x-11051015REAL PARTn202530real0-11051015IMAG PARTn202530imag0-10510 b、再取一次z=0.95e^j3pi/11，计算其DTFT，并对w绘出幅度。注意作为w的函数，幅度响应的尖峰位于何处。把尖峰位置与z的极坐标形式联系起来。 实验程序清单： b=1;

a=[1,-0.95*exp(j*3*pi/11)];

w=[-(8/11)*pi:0.05*pi:(14/11)*pi]; H=freqz(b,a,w);

stem(w,abs(H)); %以w为横轴，H的幅度为纵轴绘制离散图形

15n202530

grid;

xlabel('w'),ylabel('|H(w)|');

实验仿真图：

2018161412|H(w)|1086420-3-2-10w1234 c、如果z的角度变为3π/5，粗略绘出预期的DTFT，并对w绘出幅度，并通过对freqz的计算进行绘图来验证。

实验程序清单：

w=[-0.4*pi:0.05*pi:1.6*pi]; %W的表达式 z0=0.95*exp(j*3*pi/5);%Z0的表达式 H=1./(1-z0*exp(-j.*w));%H的表达式 subplot(211); stem(w,abs(H)); grid;

xlabel('w'),ylabel('|H(w)|'); b=1;

a=[1,-0.95*exp(j*3*pi/5)]; c=[-0.4*pi:0.05*pi:1.6*pi];

H=freqz(b,a,c); %调用freqz函数 subplot(212); %绘制子图 stem(c,abs(H)); %绘制离散图形 grid;

xlabel('c'),ylabel('|H(c)|');

实验仿真图：

2015|H(w)|1050-22015-1012w3456|H(c)|1050-2-1012c3456 d、试建立一个联系带宽与r的简单表达式。 实验程序清单：

a1=[1,-0.975*exp(j*3*pi/5)]; a2=[1,-0.95*exp(j*3*pi/5)]; a3=[1,-0.9*exp(j*3*pi/5)]; a4=[1,-0.8*exp(j*3*pi/5)]; b=1;

w=[-0.4*pi:0.1*pi:1.6*pi];

H1=freqz(b,a1,w),H2=freqz(b,a2,w),H3=freqz(b,a3,w),H4=freqz(b,a4,w); subplot(221);

stem(w,20*log(abs(H1))); grid,title('r=0.975');

xlabel('w'),ylabel('20log|H1(w)|/dB'); subplot(222);

stem(w,20*log(abs(H2))); grid,title('r=0.95');

xlabel('w'),ylabel('20log|H2(w)|/dB'); figure; subplot(223);

stem(w,20*log(abs(H3))); grid,title('r=0.9');

xlabel('w'),ylabel('20log|H3(w)|/dB'); subplot(224);

stem(w,20*log(abs(H4))); grid,title('r=0.8');

xlabel('w'),ylabel('20log|H4(w)|/dB');

实验仿真图：

r=0.96020log|H3(w)|/dBr=0.84020log|H4(w)|/dB40200-20-220002w46-20-202w46 r=0.97580606040200-20-202w46r=0.9520log|H1(w)|/dB20log|H2(w)|/dB40200-20-202w46

三：实验体会

在这次实验中，我发现了一些软件使用上的问题，比如程序运行了但是无

法显示图片等等，但是这些都在和同学交流以后得到了解决。事实上，学习是一个过程，慢慢地学到的东西肯定会更多。

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