试验离散系统的Z域分析

1 实验三、 离散系统的Z域分析 （一）实验要求 1）学习和掌握离散系统的频率特性及其幅度特性、相位特性的物理意义； 2）深入理解离散系统频率特性的对称性和周期性； 3）认识离散系统频率特性与系统参数之间的关系； 4）通过阅读、修改并调试本实验系统所给源程序，加强计算机编程能力； （二）实验内容 1、计算差分方程 （1）用MATLAB计算差分方程

当输入序列为 时的输出结果 。 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; h=filter(a,b,x); stem(k,h) xlabel('n');ylabel('h(n)') 请给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 (说明：y=filter(a,b,x)，计算系统对输入信号向量x的零状态响应输出信号向量y,x与y长度相等，其中a和b是MiiNiiinxbinya)()(所给差分方程的相量。详见教材P25-27)

2、用MATLAB计算差分方程

所对应的系统函数的FT。 差分方程所对应的系统函数为： 2

1231230.80.440.360.02()10.70.450.6zzzHzzzz





其FT为 23230.80.440.360.02()10.70.450.6jjjjjjjeeeHeeee







用MATLAB计算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度') subplot(2,2,2); plot(w/pi,imag(h));grid title('虚部') xlabel('\omega/\pi');ylabel('Amplitude') subplot(2,2,3); plot(w/pi,abs(h));grid title('幅度谱') xlabel('\omega/\pi');ylabel('幅值') subplot(2,2,4); plot(w/pi,angle(h));grid title('相位谱') xlabel('\omega/\pi');ylabel('弧度') (说明：freqz为计算数字滤波器H(z)的频率响应函数。h=freqz(num,den,w)为计算由向量w指定的数字频率点上数字滤波器H(z)的频率响应)(jeH，结果存于h向量中。Num和den为H(z)分子和分母多项式向量。详见教材P65) *练习:

①、P43中的例3、例5。 3

②、用MATLAB编程，画出,2.011)(11zzzH的频率特性图。 3、求解 6.0)6.01()5.01)(2.01(1)(21211zzzzzX

的Z反变换。 参考程序： b=1; a=poly([-0.2 0.5 0.5 -0.6 -0.6]); [r,p,k]=residuez(b,a) (说明：例

程序：b = [-4 8]; a = [1 6 8]; [r,p,k] = residuez(b,a) 运行结果：r = -12 8 p = -4 -2 k = [] 则表示：

那么： ) （三）实验报告要求 4

1、简述实验目的和实验原理，用几何确定法分析实验中选定的系统的频率特性，并与计算机计算结果相对照，根据实验结果，对系统频率特性进行讨论和总结。 2、用MATLAB编程，画出P66中的例2.6.3的频率特性图。

3、根据MATLAB求解6.0)6.01()5.01)(2.01(1)(21211zzzzzX的结

果，写出序列)(nx的表达式。 实验四、离散傅里叶变换及其快速算法 （一）实验要求 1）通过离散傅立叶变换（即DFT）的报表表示进一步了解其计算方法及意义； 2）掌握实数序列的DFT系数的对称特点； 3）学习利用DFT计算程序计算IDFT的方法； 4)学习时间抽选奇偶分解FFT算法； 5)深入掌握时间抽选奇偶分解FFT程序的编制方法；

(二)实验内容 1、对连续的单一频率周期信号 按采样频率 采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。

此时离散序列 ，即k=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下： k=8; n1=[0:1:19]; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel('t/T');ylabel('x(n)'); 5

xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel('k');ylabel('X(k)'); n2=[0:1:15]; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel('t/T');ylabel('x(n)'); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel('k');ylabel('X(k)');

6

计算结果示于图2.1，(a)和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。

2、 对如下各序列进行谱分析，绘制出其幅频特性曲线。 （1）)()(41nRnx

（2）)8()4()(2nconnconnx 对)()(41nRnx取64点FFT。由于)8()4()(2nconnconnx为周期序列，周期为16，所以取周期为16。参考程序如下： %用FFT对序列进行谱分析

%)()(41nRnx；)8()4()(2nconnconnx %X1:存放)(1nx的向量，Y1: 存放)(1kX的向量 %X2:存放)(2nx的向量，Y2: 存放)(2kX的向量 x1=[1,1,1,1,0,0,0,0]; n=0:15; x2=cos(pi*n/4)+cos(pi*n/8); i=0:7; subplot(3,2,1);stem(i,x1,'.'); axis([0 7 0 1]);%规定x轴和y轴的标值范围 xlabel('n');ylabel('x1(n)'); y1=fft(x1,8);% )(1nx 的8点FFT subplot(3,2,3);stem(i,abs(y1),'.'); xlabel('(N=8 wk=2pik/N)k');ylabel('[X1(k)]'); y1=fft(x1,64);% )(1nx 的64点FFT i=0:63; subplot(3,2,5);stem(i,abs(y1),'.'); axis([0 63 0 4]);%规定x轴和y轴的标值范围 xlabel('(N=8 wk=2pik/N)k');ylabel('[X1(k)]'); y2=fft(x2); %)(2nx的16点FFT figure;%另一幅图 subplot(2,2,1);stem(n,x2); 7

title('x2(n)的时域序列'); xlabel('n');ylabel('x2(n)'); subplot(2,2,3);stem(n,abs(y2)); title('x2(n)的幅频特性'); xlabel('(N=16 wk=2pik/N)k');ylabel('[X2(k)]');

（三）实验报告要求 1、设)]([)(),()(4nxFTeXnRnxj。分别计算)(jeX在频率区间 [0,2]上

的16点和32点等间隔采样，并绘制)(jeX采样的幅频特性图和相频特性图。 如：%DFT的MATLAB计算 xn=[1 1 1 1]; %输入时域序列向量 Xk16=fft(xn,16); %计算xn16点DFT

Xk32=fft(xn,32); %计算xn32点DFT 绘图程序略。 2、编写序列)8()4()(nconnconnx的DFT运算程序。 2解：参考程序如下： N=16; n=0:1:N-1;%时域采样 xn=cos(n*pi/4)+ cos(n*pi/8); k=0:1:N-1; %频域采样 WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^nk; %点乘方 Xk=xn*WNnk; Subplot(2,1,1) Stem(n,xn) title('xn的时域序列'); xlabel('n');ylabel('x (n)'); Subplot(2,1,2) Stem(k,abs(Xk)); title('xn的幅频特性'); xlabel('(N=16 wk=2pik/N)k');ylabel('[X (k)]');

（Fft与Ifft的说明：y=fft(x)是利用fft函数求解x的离散傅里叶变换；y=fft(x，N)，N表示离散傅里叶变换x的数据长度；函数Ifft的参数与函数