南邮dsp上机实验报告

南京邮电大学 实 验 报 告 实验名称：离散时间信号与系统的时、频域表示 离散傅立叶变换和z变换 数字滤波器的频域分析和实现 数字滤波器的设计 课程名称 数字信号处理A(双语) 班级学号________ 姓 名_____________ 开课时间 201 /201 学年， 第 学期 实验一：离散时间信号与系统的时、频域表示 一、实验目的和任务：
熟悉Matlab基本命令，理解和掌握离散时间信号与系统的时、频域表示及简单应用。在Matlab环境中，按照要求产生序列，对序列进行基本运算；
对简单离散时间系统进行仿真，计算线性时不变（LTI）系统的冲激响应和卷积输出；
计算和观察序列的离散时间傅立叶变换（DTFT）幅度谱和相位谱。

二、实验内容：
基本序列产生和运算：
Q1.1～1.3，Q1.23，Q1.30～1.33 离散时间系统仿真：
Q2.1～2.3 LTI系统：Q2.19，Q2.21，Q2.28 DTFT：Q3.1，Q3.2，Q3.4 三、实验过程与结果分析：
Q1.1 运行P1_1产生单位样本序列 u[n] 的程序与显示的波形如下: clf;n=-10:20; u=[zeros(1,10) 1 zeros(1,20)]; stem(n,u); xlabel(＇时间序号 n＇);ylabel(＇振幅＇); title(＇单位样本序列＇); axis([-10 20 0 1.2]); Q1.2 clf 命令的作用是- 清除图形窗口上的图形 axis命令的作用是- 设置坐标轴的范围和显示方式 title命令的作用是- 给图形加名字 xlabel命令的作用是- 添加x坐标标注 ylabel命令的作用是- 添加y坐标标注 Q1.3 产生有延时11个样本ud[n]的程序及其运行结果如下：
clf; n=-10:20; ud=[zeros(1,21) 1 zeros(1,9)]; stem(n,ud); xlabel(＇时间序号 n＇);ylabel(＇振幅＇); title(＇单位样本序列＇); axis([-10 20 0 1.2]); Q1.23 n=0:50; f=0.08; phase=pi/2; A=2.5; arg=2*pi*f*n-phase; x=A*cos(arg); clf; stem(n,x); axis([0 50 -3 3]); grid; title(＇正弦序列＇); xlabel(＇时间序号 n＇); ylabel(＇振幅＇); axis; 该序列的周期为12.5s Q 1.30未污染的信号s[n]是什么样的形式？加信噪声d[n]是什么的形式？ Q1.31使用语句X=s+d能产生被噪声污染的信号吗？若不能，为什么? Q1.32信号x1，x2和x3预先好x之间的关系是什么? Q1.33legend命令的作用是什么？ Q2.1 clf; n=0:100; s1=cos(2*pi*0.05*n); s2=cos(2*pi*0.47*n); x=s1+s2; M=input(＇滤波器所需的长度=＇); num=ones(1,M); y=filter(num,1,x)/M; subplot(2,2,1); plot(n,s1); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇信号#1＇); subplot(2,2,2); plot(n,s2); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇信号#2＇); subplot(2,2,3); plot(n,x); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇输入信号＇); subplot(2,2,4); plot(n,y); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇输出信号＇); axis; 输入X[n]的s[n]分量被该离散时间系统抑制。

Q2.2 程序P2.1 中 LTI system 被修改为 y[n] = 0.5(x[n]–x[n–1])后， 输入x[n] = s1[n]+s2[n] 导致的输出为: 对于输入的影响是：该系统现在是一个高通滤波器。它通过高频率的输入分量S2，而不是低频分量输入S1。

Q2.3 clf; n=0:100; s1=cos(2*pi*0.04*n); s2=cos(2*pi*0.6*n); x=s1+s2; M=input(＇滤波器所需的长度=＇); num=ones(1,M); y=filter(num,1,x)/M; subplot(2,2,1); plot(n,s1); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇信号#1＇); subplot(2,2,2); plot(n,s2); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇信号#2＇); subplot(2,2,3); plot(n,x); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇输入信号＇); subplot(2,2,4); plot(n,y); axis([0,100,-2,2]); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇输出信号＇); axis; Q2.19 运行 P2_5 生成的结果如下：: Q2.21生成的MATLAB代码如下: clf; N = 40; num = [0.9 -0.45 0.35 0.002]; den = [1.0 0.71 -0.46 -0.62]; x = [1 zeros(1,N-1)]; y = filter(num,den,x); stem(y); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); title(＇Impulse Response＇); grid; 程序产生的40个样本如下所示: Q2.28 clf; h=[3 2 1 -2 1 0 -4 0 3]; x=[1 -2 3 -4 3 2 1]; y=conv(h,x); n=0:14; subplot(2,1,1); stem(n,y); xlabel(＇时间序号n＇);ylabel(＇振幅＇); title(＇用卷积得到的输出＇);grid; x1=[x zeros(1,8)]; y1=filter(h,1,x1); subplot(2,1,2); stem(n,y1); xlabel(＇时间序号n＇);ylabel(＇振幅＇);title(＇由滤波生成的输出＇); grid; ① y[n] 和 y1[n] 的差别为：
它们无差别。

②将x[n]补零后得到 x1[n]作为输入,产生y1[n]的原因是--对于长度N1和N2的两个序列，转化率返回得到的序列长度N1 + N2-1。与此相反，过滤器接受一个输入信号和一个系统规范。返回的结果是相同的长度作为输入信号。因此，为了从转化率和滤波器得到直接比较的结果，有必要供应滤波器的输入已经零填充为长度L(x)+L(h)-1。

Q3.1 程序P3_1计算离散时间傅里叶变换的原始序列为：H(e)=; pause 命令的作用为：加参数,直接用pause的话,就是程序暂停,直至用户按任意一个按键。如果加参数，比如pause(1.5)就是程序暂停1.5秒。

Q3.2 程序 P3_1 运行结果为: DTFT 是关于的周期函数么？答：DTFT是关于的周期函数;期是 2. 四个图形的对称性为：实部是2周期和偶对称；
是2周期和奇对称；
幅度是2周期和偶对称；
相位是2周期和奇对称性。

Q3.4 修改程序 P3_1 重做Q3.2的程序如下: clf; w = -4*pi:8*pi/511:4*pi; num = [1 3 5 7 9 11 13 15 17]; den = 1; h = freqz(num, den, w); subplot(2,1,1) plot(w/pi,real(h));grid title(＇Real part of H(e^{j＼omega})＇) xlabel(＇＼omega /＼pi＇); ylabel(＇Amplitude＇); subplot(2,1,2) plot(w/pi,imag(h));grid title(＇Imaginary part of H(e^{j＼omega})＇) xlabel(＇＼omega /＼pi＇); ylabel(＇Amplitude＇); pause subplot(2,1,1) plot(w/pi,abs(h));grid title(＇Magnitude Spectrum |H(e^{j＼omega})|＇) xlabel(＇＼omega /＼pi＇); ylabel(＇Amplitude＇); subplot(2,1,2) plot(w/pi,angle(h));grid title(＇Phase Spectrum arg[H(e^{j＼omega})]＇) xlabel(＇＼omega /＼pi＇); ylabel(＇Phase in radians＇); 修改程序后的运行结果为: DTFT 是关于的周期函数么？ 答：DTFT 是关于的周期函数。周期是 - 2。

相位谱中跳变的原因是：角度返回arctan的本值。

实验名称：离散傅立叶变换和z变换 实验目的和任务：
掌握离散傅立叶变换（DFT）及逆变换（IDFT）、z变换及逆变换的计算和分析。利用Matlab语言，完成DFT和IDFT的计算及常用性质的验证，用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。

实验内容：
DFT和IDFT计算：
Q3.23～3.24 DFT的性质：
Q3.26~3.29，Q3.36，Q3.38，Q3.40 z变换分析：Q3.46～3.48 逆z变换：Q3.50 实验过程与结果分析：
Q3.23 编写一个MATLAB程序，计算并画出长度为N的L点离散傅里叶变换X[k]的值，其中L≥N，然后计算并画出L点离散傅里叶变换X[k]。对不同长度N和不同的离散傅里叶变换长度L，运行程序。讨论你的结果。

clf; N=200; L=256; nn = [0:N-1]; kk = [0:L-1]; xR = [0.1*(1:100) zeros(1,N-100)]; % real part xI = [zeros(1,N)]; x = xR + i*xI; XF = fft(x,L); subplot(3,2,1);grid; plot(nn,xR);grid; title(＇Re＼{x[n]＼}＇); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); subplot(3,2,2); plot(nn,xI);grid; title(＇Im＼{x[n]＼}＇); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); subplot(3,2,3); plot(kk,real(XF));grid; title(＇Re＼{X[k]＼}＇); xlabel(＇Frequency index k＇); ylabel(＇Amplitude＇); subplot(3,2,4); plot(kk,imag(XF));grid; title(＇Im＼{X[k]＼}＇); xlabel(＇Frequency index k＇); ylabel(＇Amplitude＇); xx = ifft(XF,L); subplot(3,2,5); plot(kk,real(xx));grid; title(＇Real part of IDFT＼{X[k]＼}＇); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); subplot(3,2,6); plot(kk,imag(xx));grid; title(＇Imag part of IDFT＼{X[k]＼}＇); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); 当N=100，L=200时，所得的图如下：
Q3.24 写一个MATLAB程序，用一个N点复数离散傅里叶计算两个长度为N的实数序列的N点离散傅里叶变换，，并将结果同直接使用两个N点离散傅里叶变换得到的结果进行比较。

clf; N=256; nn = [0:N-1]; ntime = [-N/2:N/2-1]; g = (0.75).^abs(ntime); h = (-0.9).^ntime; GF = fft(g); HF = fft(h); x = g + i*h; XF = fft(x); XFstar = conj(XF); XFstarmod = [XFstar(1) fliplr(XFstar(2:N))]; GF2 = 0.5*(XF + XFstarmod); HF2 = -i*0.5*(XF - XFstarmod); abs(max(GF-GF2)) abs(max(HF-HF2)) figure(1);clf; subplot(2,2,1);grid; plot(nn,real(GF));grid; title(＇Two N-point DFT＇＇s＇); xlabel(＇Frequency index k＇); ylabel(＇Re＼{G[k]＼}＇); subplot(2,2,2); plot(nn,imag(GF));grid; title(＇Two N-point DFT＇＇s＇); xlabel(＇Frequency index k＇); ylabel(＇Im＼{G[k]＼}＇); subplot(2,2,3);grid; plot(nn,real(GF2));grid; title(＇Single N-point DFT＇); xlabel(＇Frequency index k＇); ylabel(＇Re＼{G[k]＼}＇); subplot(2,2,4); plot(nn,imag(GF2));grid; title(＇Single N-point DFT＇); xlabel(＇Frequency index k＇); ylabel(＇Im＼{G[k]＼}＇); figure(2);clf; subplot(2,2,1);grid; plot(nn,real(HF));grid; title(＇Two N-point DFT＇＇s＇); xlabel(＇Freq index k＇); ylabel(＇Re＼{H[k]＼}＇); subplot(2,2,2); plot(nn,imag(HF));grid; title(＇Two N-point DFT＇＇s＇); xlabel(＇Freq index k＇); ylabel(＇Im＼{H[k]＼}＇); subplot(2,2,3);grid; plot(nn,real(HF2));grid; title(＇Single N-point DFT＇); xlabel(＇Freq index k＇); ylabel(＇Re＼{H[k]＼}＇); subplot(2,2,4); plot(nn,imag(HF2));grid; title(＇Single N-point DFT＇); xlabel(＇Freq index k＇); ylabel(＇Im＼{H[k]＼}＇); Q3.26 在函数circshift中，命令rem的作用是什么？ 答：rem(x,y)是用y对x求余数函数。

Q3.27 解释函数circshift怎样实现圆周移位运算。

答：在输入序列x由M的位置开始被循环移位。如果M> 0，则circshift删除从矢量x最左边开始的M个元素和它们附加在右侧的剩余元素，以获得循环移位序列。如果如果M<0，则circshift首先通过x的长度来弥补M，即序列x最右边的长度的M样品从x中删除和所附在其余的M个样本的右侧，以获得循环移位序列。

Q3.28 在函数circshift中，运算符~=的作用是什么？ 答：~=是不等于的意思。

Q3.29 解释函数circonv怎样实现圆周卷积运算。

答：输入是两个长度都为L的向量x1和x2，它是非常有用的定期延长X2的函数。让x2p成为x2延长无限长的周期的序列。从概念上讲，在定点时间上通过时序交换后的x2p的长度L交换x2p序列和x2tr等于1的元素。然后元素1至L的输出向量y是通过取x1和获得的长度为L的sh矢量之间的内积得到通过循环右移的时间反转向量x2tr。对于输出样本Y[n]的1≤N≤L时，右循环移位的量为n-1个位置上。

Q3.36 运行程序P3.9并验证离散傅里叶变换的圆周卷积性质。

g1=[1 2 3 4 5 6];g2=[1 -2 3 3 -2 1]; ycir=circonv(g1,g2); disp(＇圆周卷积的结果=＇);disp(ycir) G1=fft(g1);G2=fft(g2); yc=real(ifft(G1.*G2)); disp(＇离散傅里叶变换乘积的离散傅里叶逆变换的结果=＇);disp(yc) function y=circonv(x1,x2) L1=length(x1);L2=length(x2); if L1~=L2,error(＇长度不相等的序列＇),end y=zeros(1,L1); x2tr=[x2(1) x2(L2:-1:2)]; for k=1:L1 sh=circshift(x2tr,1-k); h=x1.*sh; y(k)=sum(h); end 圆周卷积的结果= 12 12 12 12 12 12 离散傅里叶变换乘积的离散傅里叶逆变换的结果= 12 28 14 0 16 14 Q3.38 运行程序P3.10并验证线性卷积可通过圆周卷积得到。

Linear convolution via circular convolution = 2 6 10 15 21 15 7 9 5 Direct linear convolution = 2 6 10 15 21 15 7 9 5 使用圆周卷积确实有可能得到线性卷积 Q3.40 编写一个MATLAB程序，对两个序列做离散傅里叶变换，已生成他们的线性卷积。用此程序验证Q3.38和Q3.39的结果 编写的MATLAB程序：
g1 = [1 2 3 4 5]; g2 = [2 2 0 1 1]; g1e = [g1 zeros(1,length(g2)-1)]; g2e = [g2 zeros(1,length(g1)-1)]; G1EF = fft(g1e); G2EF = fft(g2e); ylin = real(ifft(G1EF.*G2EF)); disp(＇Linear convolution via DFT = ＇); disp(ylin); Linear convolution via DFT = 2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.0000 9.0000 5.0000 Q3.46 使用程序P3.1在单位圆上求下面的z变换：
G(z)= 使用函数zplane可以很容易地得到有理Z变换G(z)的极零点图，该函数有两种形式，若Z变换用式（3.32）所示的有理函数的形式给出，使用的命令是zplane(num,den), 其中num和den是按z^-1的升幂排列的G（z）的分子和分母多项式系数的行向量，另外，若给出G（z）的零点和极点，将用到的命令是zplane（zeros，poles），其中zeros和poles都是列向量，在由matlab产生的极零点图中，几点的位置用符号x表示，而零点的位置符号用o表示。

得到结果如图：
Q3.47 编写一个MATLAB程序，计算并显示零点和极点，计算并显示其因式形式，并产生z的两个多项式之比的形式表示的z变换的极零点图。使用该程序，分析式(3.32)的z变换G(z)。

编写的MATLAB程序：
clf; % initialize num = [2 5 9 5 3]; den = [5 45 2 1 1]; [z p k] = tf2zpk(num,den); disp(＇Zeros:＇); disp(z); disp(＇Poles:＇); disp(p); input(＇Hit <return> to continue...＇); [sos k] = zp2sos(z,p,k) input(＇Hit <return> to continue...＇); zplane(z,p); 运行结果：
Zeros: -1.0000 + 1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614i Poles: -8.9576 -0.2718 0.1147 + 0.2627i 0.1147 - 0.2627i sos = 1.0000 2.0000 3.0000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822 k = 0.4000 Q3.48 通过习题Q3.47产生的极零点图，求出G(z)的收敛域的数目。清楚地显示所有的收敛域。由极零点图说明离散时间傅里叶变换是否存在。

R1 : | z | < 0.2718 (left-sided, not stable) R2 : 0.2718 < | z | < 0.2866 (two-sided, not stable) R 3: 0.2866 < | z | < 8.9576 (two-sided, stable) R4 : | z | > 8.9576 (right-sided, not stable) 不能从极零点图肯定的说DTFT是否存在，因为其收敛域一定要指定。当收敛域在上述R 3内所获得的序列却是证明了DTFT的存在，它是一个具有双面冲激响应的稳定系统。

Q3.50 编写一个MATLAB程序，计算一个有理逆z变换的前L个样本，其中L的值由用户通过命令input提供。用该程序计算并画出式(3.32)中G(z)的逆变换的前50个样本。使用命令stem画出由逆变换产生的序列。

编写的MATLAB程序：
clf; num = [2 5 9 5 3]; den = [5 45 2 1 1]; L = input(＇Enter the number of samples L: ＇); [g t] = impz(num,den,L); stem(t,g); title([＇First ＇,num2str(L),＇ samples of impulse response＇]); xlabel(＇Time Index n＇); ylabel(＇h[n]＇); Enter the number of samples L: 50 实验名称：数字滤波器的频域分析和实现 一、实验目的和任务：
（1）求滤波器的幅度响应和相位响应，观察对称性，判断滤波器类型。

（2）用Matlab函数实现系统的级联型和并联型结构，并对滤波器进行结构仿真。

二、实验内容：
系统传递函数的级联和并联实现：Q6.1, Q6.3, Q6.5, Q8.3,Q8.5, 滤波器的幅频特性分析：Q8.1, Q8.9,Q8.10,Q8.14 三、实验过程与结果分析：
Q6.1 使用程序P6.1，生成如下有限冲激响应传输函数的一个级联实现：
H(z)=2+10z+23z+34z+31z+16z 编写的MATLAB程序：
% Program P6_1 num = input(＇Numerator coefficient vector = ＇); den = input(＇Denominator coefficient vector = ＇); [b,a] = eqtflength(num,den); % make lengths equal [z,p,k] = tf2zp(b,a); sos = zp2sos(z,p,k) Numerator coefficient vector = [2 10 23 34 31 16 4] Denominator coefficient vector = [1] sos = 2.0000 6.0000 4.0000 1.0000 0 0 1.0000 1.0000 2.0000 1.0000 0 0 1.0000 1.0000 0.5000 1.0000 0 0 画出级联实现的框图：
H(z)不是一个线性相位传输函数。

Q6.3 使用程序P6.1，H(z)= Numerator coefficient vector = [3 8 12 7 2 –2] Denominator coefficient vector = [16 24 24 14 5 1] sos = 0.1875 -0.0625 0 1.0000 0.5000 0 1.0000 2.0000 2.0000 1.0000 0.5000 0.2500 1.0000 1.0000 1.0000 1.0000 0.5000 0.5000 画出级联实现的框图：
Q6.5 使用程序P6.2，式H(z)=所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图：
Parallel Form I Residues are -0.4219 + 0.6201i -0.4219 - 0.6201i 2.3437 0.3437 - 2.5079i 0.3437 + 2.5079i Poles are at -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i Constant value -2 Parallel Form II Residues are -0.3047 - 0.4341i -0.3047 + 0.4341i -1.1719 1.0000 + 0.7758i 1.0000 - 0.7758i Poles are at -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i Constant value 0.1875 Q8.1 程序P8.1设计的滤波器类型为 IIR带阻滤波器 ，其指标为 =[0.4 0.5] =[0.1 0.8] , R=1dB，R=30dB，滤波器的阶数为 4 。为了验证仿真需计算的冲激响应样本为 5个。仿真是正确的吗？ % Program P8_1 Wp = [0.4 0.5]; Ws = [0.1 0.8]; Rp = 1; Rs = 30; [N1, Wn1] = buttord(Wp, Ws, Rp, Rs); [num,den] = butter(N1,Wn1); disp(＇Numerator coefficients are ＇);disp(num); disp(＇Denominator coefficients are ＇);disp(den); impres = direct2(num,den,[1 zeros(1,4*N1)]); [p,d] = strucver(impres,2*N1); disp(＇Actual numerator coeffs are ＇); disp(p＇); disp(＇Actual denominator coeffs are ＇); disp(d＇); Numerator coefficients are 0.0571 0 -0.1143 0 0.0571 Denominator coefficients are 1.0000 -0.5099 1.2862 -0.3350 0.4479 Actual numerator coeffs are 0.0571 -0.0000 -0.1143 0.0000 0.0571 Actual denominator coeffs are 1.0000 -0.5099 1.2862 -0.3350 0.4479 直接II型结构的因果无限冲激响应滤波器，滤波器的阶数是4阶 Q8.3 生成Q8.1中产生的传输函数的一个级联实现，并编写一个程序来仿真它。其中每个单独的部分用一个直接Ⅱ型实现。验证仿真。

format short num = input(＇Numerator coefficients = ＇); den = input(＇Denominator coefficients = ＇); Numfactors = factorize(num); Denfactors = factorize(den); disp(＇Numerator Factors＇),disp(Numfactors) disp(＇Denominator Factors＇),disp(Denfactors) Numerator coefficients = [0.0571 0 -0.1143 0 0.0571] Denominator coefficients = [1.0000 -0.5099 1.2862 -0.3350 0.4479] Numerator Factors Columns 1 through 3 1.000000000000000 1.021143238995107 0 1.000000000000000 0.979294541463240 0 1.000000000000000 -1.021143238995070 0 1.000000000000000 -0.979294541463277 0 Denominator Factors Columns 1 through 3 1.000000000000000 -0.597633166080704 0.678485370310279 1.000000000000000 0.087733166080704 0.660146879504817 Q8.5 生成Q8.1中传输函数的一个并联Ⅰ型实现，并编写一个程序来仿真它。每一个单独的部分用直接Ⅱ型实现。验证仿真。

num = input(＇Numerator coefficient vector = ＇); den = input(＇Denominator coefficient vector = ＇); [r1,p1,k1] = residuez(num,den); [r2,p2,k2] = residue(num,den); disp(＇Parallel Form I＇) disp(＇Residues are＇);disp(r1); disp(＇Poles are at＇);disp(p1); disp(＇Constant value＇);disp(k1); disp(＇Parallel Form II＇) disp(＇Residues are＇);disp(r2); disp(＇Poles are at＇);disp(p2); disp(＇Constant value＇);disp(k2); [b1,a1]=residuez(R1,P1,0); R1=[r1(1) r1(2) r1(3) r1(4)]; P1=[p1(1) p1(2) p1(3) p1(4)]; disp(＇b1=＇);disp(b1); disp(＇a1=＇);disp(a1); Numerator coefficient vector = [0.0571 0 -0.1143 0 0.0571] Denominator coefficient vector = [1.0000 -0.5099 1.2862 -0.3350 0.4479 ] Parallel Form I Residues are -0.023475623417519 + 0.197736845648730i -0.023475623417519 - 0.197736845648730i -0.011716283258078 - 0.213037612139943i -0.011716283258078 + 0.213037612139943i Poles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311i Constant value 0.127483813351194 Parallel Form II Residues are -0.158795680938490 + 0.041067400713457i -0.158795680938490 - 0.041067400713457i 0.173353325938490 - 0.000160296180734i 0.173353325938490 + 0.000160296180734i Poles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311i Constant value 0.057100000000000 b1= -0.046951246835039 -0.289531739579519 0 a1= 1.000000000000000 -0.597633166080704 0.678485370310279 Q8.9 程序P8.3设计的滤波器类型为 IIR低通椭圆滤波器 ，指标为=0.25 =0.55，R=0.5dB，R=50dB，阶数为 8，输入的正弦序列的频率为f=0.7kHz。

% Program P8_3 % Illustration of Filtering by an IIR Filter % clf; % Generate the input sequence k = 0:50; w2 = 0.7*pi;w1 = 0.2*pi; x1 = 1.5*cos(w1*k); x2 = 2*cos(w2*k); x = x1+x2; % Determine the filter transfer function [N, Wn] = ellipord(0.25, 0.55, 0.5, 50); [num, den] = ellip(N,0.5, 50,Wn); % Generate the output sequence y = filter(num,den,x); % Plot the input and the output sequences subplot(2,1,1); stem(k,x); grid; axis([0 50 -4 4]); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); title(＇Input Sequence＇); subplot(2,1,2); stem(k,y); grid; axis([0 50 -4 4]); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); title(＇Output Sequence＇); Q8.10 运行P8.3并产生两个图形。哪种输入成分会在滤波器输出出现？为什么输出序列的开始部分不是一种理想的正弦曲线？修改P8.3，以便纸过滤序列X2[n]。产生的输出序列和预料的是一样的吗？ 答：产生的输出序列和预料的是不一样的 Q8.14 程序P8.4设计的滤波器类型为 FIR低通滤波器，指标为=[0 0.3] =[0.5 ] ，阶数为 9 ，为了验证仿真需计算的冲激响应样本为 10个。仿真是正确的吗？ 答：滤波器类型为 FIR低通滤波器，指标为=[0 0.3] =[0.5 ] ，阶数为 9 ，为了验证仿真需计算的冲激响应样本为 10个，仿真是正确的。

实验名称：数字滤波器的设计 一、实验目的和任务：
（1）用窗口法设计满足指标的FIR数字滤波器。

（2）以模拟低通滤波器为原型设计IIR数字滤波器 （3）选定一个信号滤波问题，设计数字滤波器，验证滤波效果。

二、实验内容：
阅读Page 91-93相应的函数和程序P7.1,完成Q7.1, Q7.5,Q7.6 阅读Page 94-96相应的函数, 完成Q7.9,Q7.13,Q7.14,Q7.20 sinc函数的功能与使用，可通过help查询：
Matlab-Help-Search-Function Name-输入sinc 实验指导书Page 49有sinc函数使用实例 幅度响应的分析：
通过DTFT定义计算幅度响应 通过freqz函数分析 三、实验过程与结果分析：
Q7.1 用MATLAB确定一个数字无线冲激响应低通滤波器所有四种类型的最低阶数。指标如下：40kHz的抽样率，4kHZ的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。评论你的结果。

根据题意：F=40kHz，F=4 kHz，F=8kHz，通带增益R=0.5dB，阻带增益R=40dB 可以得出= W = W (1) 根据上述数据和buttord函数[N, Wn] = buttord(0.2,0.4,0.5,40)得到巴特沃兹滤波器的最低阶数N=8。

Wn=0.2469 =0.2469 (2) 根据上述数据和cheb1ord函数[N, Wn] = cheb1ord(0.2,0.4,0.5,40)得到切比雪夫1型滤波器的最低阶数N=5。

Wn=0.2000 =0.2000 (3) 根据上述数据和cheb2ord函数[N, Wn] = cheb1ord(0.2,0.4,0.5,40)得到切比雪夫2型滤波器的最低阶数N=5。

Wn=0.4000 =0.4000 (4) 根据上述数据和ellipord函数[N, Wn] = ellipord (0.2,0.4, 0.5,40)得到椭圆滤波器的最低阶数N=4。

Wn=0.2000 =0.2000 Q7.5 通过运行P7.1设计巴特沃兹带阻滤波器，写出所产生的传输函数的准确表达式。滤波器的指标是什么？使用MATLAB计算并绘制滤波器未畸变的相位相应及群延迟相应。

MATLAB程序为：
% Program P7_1 % Design of a Butterworth Bandstop Digital Filter Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50; % Estimate the Filter Order [N1, Wn1] = buttord(0.2, 0.4, 0.5, 40); % Design the Filter [num,den] = butter(N1,Wn1,＇stop＇); % Display the transfer function disp(＇Numerator Coefficients are ＇);disp(num); disp(＇Denominator Coefficients are ＇);disp(den); % Compute the gain response [g, w] = gain(num,den); % Plot the gain response plot(w/pi,g);grid axis([0 1 -60 5]); xlabel(＇＼omega /＼pi＇); ylabel(＇Gain in dB＇); title(＇Gain Response of a Butterworth Bandstop Filter＇); 得到：
Numerator Coefficients are Columns 1 through 9 0.0493 0.0000 0.2465 0.0000 0.4930 0.0000 0.4930 0.0000 0.2465 Columns 10 through 11 0.0 0.0493 Denominator Coefficients are Columns 1 through 9 1.0000 0.0000 -0.0850 0.0000 0.6360 0.0000 -0.0288 0.0000 0.0561 Columns 10 through 11 0.0000 -0.0008 因此表达式为：
H (z)= 指标为=0.2，=0.4，=0.6，=0.8，R=0.4dB，R=50dB Q7.6 修改P7.1后来设计符合习题Q7.1所给指标的切比雪夫1型低通滤波器。写出产生的传输函数的准确表达式。使用MATLAB计算并绘制滤波器未畸变的相位相应及群延迟响应。

修改P7.1后的程序（切比雪夫1型低通滤波器）：
% Program Q7_6 % Design spec as given in Q7.1. FT = 40*10^3; % sampling freq Fp = 4*10^3; % analog passband edge freq Fs = 8*10^3; % analog stopband edge freq Rp = 0.5; % max passband ripple, dB Rs = 40; % min stopband attenuation, dB % Convert spec to normalized digital frequencies omega_p = 2*pi*Fp/FT; Wp = 2*Fp/FT; % omega_p/pi omega_s = 2*pi*Fs/FT; Ws = 2*Fs/FT; % omega_s/pi % Estimate the Filter Order [N, Wn] = cheb1ord(Wp, Ws, Rp, Rs); % Design the Filter [num,den] = cheby1(N,Rp,Wn); % Display the transfer function disp(＇Numerator Coefficients are ＇);disp(num); disp(＇Denominator Coefficients are ＇);disp(den); % Compute the gain response [g, w] = gain(num,den); % Plot the gain response figure(1); plot(w/pi,g);grid; axis([0 1 -60 5]); xlabel(＇＼omega /＼pi＇); ylabel(＇Gain in dB＇); title(＇Gain Response of a Type 1 Chebyshev Lowpass Filter＇); % Find and plot the phase figure(2); w2 = 0:pi/511:pi; Hz = freqz(num,den,w2); Phase = unwrap(angle(Hz)); plot(w2/pi,Phase);grid; xlabel(＇＼omega /＼pi＇); ylabel(＇Unwrapped Phase (rad)＇); title(＇Unwrapped Phase Response of a Type 1 Chebyshev Lowpass Filter＇); % Find and plot the group delay figure(3); GR = grpdelay(num,den,w2); plot(w2/pi,GR);grid; xlabel(＇＼omega /＼pi＇); ylabel(＇Group Delay (sec)＇); title(＇Group Delay of a Type 1 Chebyshev Lowpass Filter＇); 编写的MATLAB程序如下（计算未畸变的相位响应及群延迟响应）：
% Program Q7_6 Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50; % Estimate the Filter Order [N1, Wn1] = buttord(Wp, Ws, Rp, Rs); % Design the Filter [num,den] = butter(N1,Wn1,＇stop＇); % Find the frequency response; find and plot unwrapped phase wp = 0:pi/1023:pi; wg = 0:pi/511:pi; Hz = freqz(num,den,wp); Phase = unwrap(angle(Hz)); figure(1); plot(wp/pi,Phase); grid; % axis([0 1 a b]); xlabel(＇＼omega /＼pi＇); ylabel(＇Unwrapped Phase (rad)＇); title(＇Unwrapped Phase Response of a Butterworth Bandstop Filter＇); % Find and plot the group delay GR = grpdelay(num,den,wg); figure(2); plot(wg/pi,GR); grid; %axis([0 1 a b]); xlabel(＇＼omega /＼pi＇); ylabel(＇Group Delay (sec)＇); title(＇Group Delay of a Butterworth Bandstop Filter＇); Q7.9 使用sinc编写一个MATLAB程序，以产生截止频率在=0.4处，长度分别为81,61,41,和21的四个零相位低通滤波器的冲激响应系数，然后计算并画出它们的幅度响应。使用冒号“：”运算符从长度为81的滤波器的冲激响应系数中抽出较短长度滤波器的冲激响应系数。在每一个滤波器的截止频率两边研究频率响应的摆动行为。波纹的数量与滤波器的长度有什么关系？最大波纹的高度与滤波器的长度有什么关系？将怎样修改上述程序以产生一个偶数长度的零相位低通滤波器的冲激响应系数？ （1）长度为81的零相位低通滤波器的冲激响应：
（2）长度为61的零相位低通滤波器的冲激响应：
（3）长度为41的零相位低通滤波器的冲激响应：
（4）长度为21的零相位低通滤波器的冲激响应：
答：从这些图形可以证明不同情况下的幅值振动都是由于Gibb效应。波纹的数量的减少和滤波器的长度成正比的。最大波纹的高度与滤波器的长度没有关系。

修改程序产生一个偶数长度的零相位低通滤波器的冲激响应系数：
% Program Q7_9 n = -39.5:39.5; % this gives us a length of 80 hn_80 = 0.4 * sinc(0.4*n); % the length-80 impulse response omega = 0:pi/1023:pi; % radian frequency vector W = omega/pi; % Matlab normalized freq vector Hz_80 = abs(freqz(hn_80,1,omega)); % 1024 samles of |H(e^jw)| figure(1); plot(W,Hz_80); grid; xlabel(＇＼omega /＼pi＇); ylabel(＇|H(e^{j＼omega})|＇); title(＇Magnitude Response for Length=80＇); % Reduce length to 60 and repeat hn_60 = hn_80(11:70); Hz_60 = abs(freqz(hn_60,1,omega)); figure(2); plot(W,Hz_60); grid; xlabel(＇＼omega /＼pi＇); ylabel(＇|H(e^{j＼omega})|＇); title(＇Magnitude Response for Length=60＇); % Reduce length to 40 and repeat hn_40 = hn_60(11:50); Hz_40 = abs(freqz(hn_40,1,omega)); figure(3); plot(W,Hz_40); grid; xlabel(＇＼omega /＼pi＇); ylabel(＇|H(e^{j＼omega})|＇); title(＇Magnitude Response for Length=40＇); % Reduce length to 20 and repeat hn_20 = hn_40(11:30); Hz_20 = abs(freqz(hn_20,1,omega)); figure(4); plot(W,Hz_20); grid; xlabel(＇＼omega /＼pi＇); ylabel(＇|H(e^{j＼omega})|＇); title(＇Magnitude Response for Length=20＇); Q7.13使用函数kaiord，估计具有以下指标的线性相位低通有限冲激相应滤波器的阶数：通带边界为2kHz，阻带边界为2.5kHz，通带波纹，阻带波纹=0.005.抽样率为10kHz。在函数kaiord中命令ceil和nargin的作用是什么？ 答：使用函数kaiord和上述数据可以估计得到N=46 命令ceil：向正方向舍入成整数的函数。

命令nargin：用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。通常可以用他来设定一些默认值。

Q7.14 对下面的情况重做Q7.13：(a)20kHz的抽样率，(b) 0.002和=0.002，(c)阻带边界为2.3kHz。把每一种情况中得到的滤波器长度与Q7.13中得到的滤波器长度相比作比较。评论抽样率、波纹以及过渡带宽对滤波器阶数的影响。

(a)20kHz的抽样率时，N=91；

(b) 0.002和=0.002时，N=57；

(c)阻带边界为2.3kHz时，N=76；

1、抽样率对滤波器阶数的影响：对于一个给定的模拟过渡带宽，抽样率的增长导致估计的滤波器阶数的成比例增长，一直到下一个整数值。

2、波纹对滤波器阶数的影响：估计的滤波器阶数大约与波纹的log值成一定比例。

3、过渡带宽对滤波器阶数的影响：在一定的范围内滤波器阶数与过渡带宽成比例的变化。

课后习题：
M 2.2 (a)Using Program2_2,generate the sequence shown in Figures 2.23 and 2.24 (b)Generate and plot the complex exponential sequence -2.7e^((-0.4+jπ/6)n) for 0≤n≤82 using Program 2_2. Answer: (a) For figure 2.23: a =-1/12; b = pi/6; c = a + b*i;  K = 1; N =41;  n = 1:N;  x = K*exp(c*n);  stem(n,real(x));   xlabel(＇时间n＇);ylabel(＇振幅＇);  title(＇实数部分＇);  disp(＇PRESS RETURN for imaginary part＇);  pause  stem(n,imag(x));%Plot the imaginary part  xlabel(＇Time index n＇);ylabel(＇Amplitude＇);  title(＇Imaginary part＇); For figure 2.24(a): a =log(1.2); b = 0; c = a + b*i;  K = 0.2; N =31;  n = 1:N;  x = K*exp(c*n);  stem(n,real(x));   xlabel(＇时间n＇);ylabel(＇振幅＇);  title(＇实数部分＇);  disp(＇PRESS RETURN for imaginary part＇);  pause  stem(n,imag(x));%Plot the imaginary part  xlabel(＇Time index n＇);ylabel(＇Amplitude＇);  title(＇Imaginary part＇); For figure 2.24(b): a =log(0.9); b = 0; c = a + b*i;  K = 20; N =31;  n = 1:N;  x = K*exp(c*n);  stem(n,real(x));   xlabel(＇时间n＇);ylabel(＇振幅＇);  title(＇实数部分＇);  disp(＇PRESS RETURN for imaginary part＇);  pause  stem(n,imag(x));%Plot the imaginary part  xlabel(＇Time index n＇);ylabel(＇Amplitude＇);  title(＇Imaginary part＇); (b) a =-0.4; b = pi/6; c = a + b*i;  K = -2.7; N =8;  n = 1:N;  x = K*exp(c*n);  stem(n,real(x));   xlabel(＇时间n＇);ylabel(＇振幅＇);  title(＇实数部分＇);  disp(＇PRESS RETURN for imaginary part＇);  pause  stem(n,imag(x));%Plot the imaginary part  xlabel(＇Time index n＇);ylabel(＇Amplitude＇);  title(＇Imaginary part＇); M2-4 (a) L=input(＇desired length = ＇) A=input(＇Amplitude = ＇) omega=input(＇angular frequency = ＇) phi = input(＇phase = ＇) n=0:L-1 x=A*cos(omega*n+phi) stem(n,x) xlabel(＇time index＇ ) ylabel(＇amplitude＇) title([＇＼omega_{o} = ＇,num2str(omega/pi),＇＼pi＇]) (b) M2.6 Write a MATLAB program to plot a continuous-time sinusoidal signal and its sampled version, and verify Figure 2.28. You need to use the hold function to keep both plots. t=0:0.001:1; fo=input(＇Frequency of sinusoid in Hz=＇); FT=input(＇Sampling frequency in Hz=＇); g1=cos(2*pi*fo*t); plot(t,g1,＇-＇); xlabel(＇time＇); ylabel(＇Amplitude＇); hold n=0:1:FT; gs=cos(2*pi*fo*n/FT); plot(n/FT,gs,＇o＇); hold off Frequency of sinusoid in Hz=10 Sampling frequency in Hz=5 clf; w=-4*pi:8*pi/511:4*pi; num=[0.1323 0.01910412 -0.05978637 0.01910412 0.1323]; den=[1 0.1386 0.8258 0.1393 0.4153]; h=freqz(num,den,w); subplot(2,1,1) plot(w/pi,real(h));grid title(＇(e^{j＼omega})的实部＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇) subplot(2,1,2) plot(w/pi,imag(h));grid title(＇H(e^{j＼omega}的虚部)＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); pause subplot(2,1,1) plot(w/pi,abs(h));grid title(＇H(e^{j＼omega})的实部＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); subplot(2,1,2) plot(w/pi,angle(h));grid title(＇[H(e^{j＼omega})]的虚部＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); clf; w=-4*pi:8*pi/511:4*pi; num=[0.3192 0.0601692 -0.0601692 -0.3192]; den=[1 0.7856 1.4654 -0.2346]; h=freqz(num,den,w); subplot(2,1,1) plot(w/pi,real(h));grid title(＇(e^{j＼omega})幅度谱＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇) subplot(2,1,2) plot(w/pi,imag(h));grid title(＇H(e^{j＼omega}相位谱)＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); pause subplot(2,1,1) plot(w/pi,abs(h));grid title(＇H(e^{j＼omega})幅度谱＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); subplot(2,1,2) plot(w/pi,angle(h));grid title(＇相位谱[H(e^{j＼omega})]＇) xlabel(＇＼omega/＼pi＇); ylabel(＇振幅＇); M6.5 Repeat Problem 6.72 using Matlab A casual stable LTI discrete-time system is characterized by an impulse response h1[n]=1.9δ[n]+0.5(-0.2)^nu[n]-0.6(0.7)^nu[n] Determine the impulse response h2[n] of its inverse system,which is causal and stable. h1= Columns 1 through 10 1.8000 -0.5200 -0.2740 -0.2098 -0.1433 -0.1010 -0.0706 -0.0496 -0.0346 -0.0242 Columns 11 through 20 -0.0169 -0.0119 -0.0083 -0.0058 -0.0041 -0.0028 -0.0020 -0.0014 -0.0010 -0.0007 The first 20 samples of h2[n] are h2= Columns 1 through 10 0.5556 0.1605 0.1310 0.1270 0.1196 0.1131 0.1069 0.1010 0.0955 0.0903 Columns 11 through 20 0.0853 0.0806 0.0762 0.0720 0.0681 0.0644 0.0608 0.0575 0.0544 0.0514 The first 20 samples of the convolution are Columns 1 through 10 1.0001 0.0001 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 Columns 11 through 20 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.000 8-2 clf; k = 0:50; w2 = 0.7*pi;w1 = 0.2*pi; x1 = 1.5*cos(w1*k); x2 = 2*cos(w2*k); x = x1+x2; [N, Wn] = ellipord(0.25, 0.55, 0.5, 50); [num, den] = ellip(N,0.5, 50,Wn); y = filter(num,den,x); subplot(2,1,1); stem(k,x); grid; axis([0 50 -4 4]); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); title(＇Input Sequence＇); subplot(2,1,2); stem(k,y); grid; axis([0 50 -4 4]); xlabel(＇Time index n＇); ylabel(＇Amplitude＇); title(＇Output Sequence＇); 2.1 2.2 3.3 9.1 数字信号处理实验小结及心得体会：