南邮2013DSP实验报告(英文mitra版)

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

2、算： Q1.11.3，Q1.23，Q1.301.33离散时间系统仿真： Q2.12.3LTI系统：Q2.19，Q2.21，Q2.28DTFT：Q3.1，Q3.2，Q3.4实验过程与结果分析：Q1.1 运行P1_1产生单位样本序列 un 的程序与显示的波形如下: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

3、命令的作用是- 设置坐标轴的范围和显示方式title命令的作用是- 给图形加名字xlabel命令的作用是- 添加x坐标标注ylabel命令的作用是- 添加y坐标标注Q1.3产生有延时11个样本udn的程序及其运行结果如下： clf;n=-10:20;u=zeros(1,21) 1 zeros(1,9);stem(n,u);xlabel('时间序号n');ylabel('振幅');title('单位样本序列');axis(-10 20 0 1.2);Q1.23 修改上述程序，以长生长度为50、频率为0.08、振幅为2.5、相移为90度的一个正弦序列

4、并显示它。该序列的周期是多少？n=0:50;f=0.1;f=0.08;phase=pi;A=2.5;arg=2*pi*f*n-phase;x=A*cos(arg);clf;stem(n,x);axis(0 50 -3 3);grid;周期为：2/ = 1/f = 1/0.08 = 1/(8/100) = 100/8 = 25/2.Q1.30 未污染的信号sn 是什么样的形式？加性噪声dn 是什么样的形式？加性噪声dn 是均匀分布在 -0.4和+0.4之间的随机序列Q1.31 使用语句s=s+d能产生被噪声污染的信号吗？若不能，为什么？不能。因为d是列向量，s是行向量Q1.32 信号x1、x2、

5、x3与x之间的关系是什么？x1是x的延时，x2和x相等，x3超前于xQ1.33 legend的作用是什么legend用于产生图例说明Q2.1 对于M = 2 和输入 xn = s1n+s2n，程序P2.1的输出为： 输入 xn 被该离散时间系统抑制的分量为- Signal #2的高频分量 Q2.2 程序P2.1 中 LTI system 被修改为 yn = 0.5(xnxn1)后， 输入 xn = s1n+s2n 导致的输出为:对于输入的影响是- 该系统现在是一个高通滤波器。它通过高频率的输入分量S2，而不是低频分量输入S1。Q2.3程序 P2_1 对于不同M（M=4,6）取值和不同正弦分量(

6、任取2个)取值的运行结果如下： M=4 f1=0.05 f2=0.10M=6 f1=0.30 f2=0.30Q2.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('时间序号n'); ylabel('振幅');title('冲击响应'); grid;程序产生的40个样本如下所示:

7、Q2.28 程序P2_7产生的序列 yn and y1n 如下所示： yn 和 y1n 的差别为- 它们无差别。将xn补零后得到 x1n作为输入,产生y1n的原因是-对于长度N1和N2的两个序列，转化率返回得到的序列长度N1 + N2-1。与此相反，过滤器接受一个输入信号和一个系统规范。返回的结果是相同的长度作为输入信号。因此，为了从转化率和滤波器得到直接比较的结果，有必要供应滤波器的输入已经零填充为长度L(x)+L(h)-1。Q3.1 程序P3_1计算离散时间傅里叶变换的原始序列为-H(e)=pause 命令的作用为- 不加参数,直接用pause的话,就是程序暂停,直至用户按任意一个按键。如

8、果加参数，比如pause(1.5)就是程序暂停1.5秒。Q3.2 程序 P3_1 运行结果为:DTFT 是关于 w的周期函数么？答：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);gridtitle(&#

9、39;H(ejomega)的实部')xlabel('omega /pi');ylabel('振幅');subplot(2,1,2)plot(w/pi,imag(h);gridtitle('H(ejomega)的虚部')xlabel('omega /pi');ylabel('振幅');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle('|H(ejomega)|幅度谱')xlabel('omega /pi');ylabel('振幅

10、');subplot(2,1,2)plot(w/pi,angle(h);gridtitle('相位谱argH(ejomega)')xlabel('omega /pi');ylabel('以弧度为单位的相位');修改程序后的运行结果为: DTFT 是关于 w的周期函数么？ 答：DTFT 是关于 w的周期函数。周期是 - 2 相位谱中跳变的原因是 - 角度返回到arctan的本值实验名称：离散傅立叶变换和z变换实验目的和任务：掌握离散傅立叶变换（DFT）及逆变换（IDFT）、z变换及逆变换的计算和分析。利用Matlab语言，完成DFT和IDF

11、T的计算及常用性质的验证，用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。实验内容：DFT和IDFT计算： Q3.233.24DFT的性质： Q3.263.29，Q3.36，Q3.38，Q3.40z变换分析：Q3.463.48逆z变换：Q3.50实验过程与结果分析：Q3.23 编写一个MATLAB程序，计算并画出长度为N的L点离散傅里叶变换Xk的值，其中LN，然后计算并画出L点离散傅里叶变换Xk。对不同长度N和不同的离散傅里叶变换长度L，运行程序。讨论你的结果。编写的MATLAB程序：clf;N=200; % length of signalL=256; % length of

12、DFTnn = 0:N-1;kk = 0:L-1;xR = 0.1*(1:100) zeros(1,N-100); xI = zeros(1,N); x = xR + i*xI;XF = fft(x,L);subplot(3,2,1);grid;plot(nn,xR);grid;title('Rexn');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,2);plot(nn,xI);grid;title('Imxn');xlabel('Time index n&

13、#39;);ylabel('Amplitude');subplot(3,2,3);plot(kk,real(XF);grid;title('ReXk');xlabel('Frequency index k');ylabel('Amplitude');subplot(3,2,4);plot(kk,imag(XF);grid;title('ImXk');xlabel('Frequency index k');ylabel('Amplitude');% IDFTxx = ifft(XF,L

14、);subplot(3,2,5);plot(kk,real(xx);grid;title('Real part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,6);plot(kk,imag(xx);grid;title('Imag part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');Q3.24 写一个MATLAB程序，用一个N点复数离散傅里叶计算

15、两个长度为N的实数序列的N点离散傅里叶变换，并将结果同直接使用两个N点离散傅里叶变换得到的结果进行比较。 编写的MATLAB程序：clf;N=256; % length of signalnn = 0:N-1;ntime = -N/2:N/2-1;g = (0.75).abs(ntime); % signal gh = (-0.9).ntime; % signal hGF = fft(g);HF = fft(h);x = g + i*h;XF = fft(x);XFstar = conj(XF);XFstarmod = XFstar(1) fliplr(XFstar(2:N);GF2 = 0.

16、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('ReGk');subplot(2,2,2);plot(nn,imag(GF);grid;title('Two N-point DF

17、T''s');xlabel('Frequency index k');ylabel('ImGk');subplot(2,2,3);grid;plot(nn,real(GF2);grid;title('Single N-point DFT');xlabel('Frequency index k');ylabel('ReGk');subplot(2,2,4);plot(nn,imag(GF2);grid;title('Single N-point DFT');xlabel(

18、9;Frequency index k');ylabel('ImGk');figure(2);clf;subplot(2,2,1);grid;plot(nn,real(HF);grid;title('Two N-point DFT''s');xlabel('Freq index k');ylabel('ReHk');subplot(2,2,2);plot(nn,imag(HF);grid;title('Two N-point DFT''s');xlabel('Freq

19、 index k');ylabel('ImHk');subplot(2,2,3);grid;plot(nn,real(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ReHk');subplot(2,2,4);plot(nn,imag(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ImHk');

20、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 解释函数circo

21、nv怎样实现圆周卷积运算。 答：输入是两个长度都为L的向量x1和x2，它是非常有用的定期延长X2的函数。让x2p成为x2延长无限长的周期的序列。从概念上讲，在定点时间上通过时序交换后的x2p的长度L交换x2p序列和x2tr等于1的元素。然后元素1至L的输出向量y是通过取x1和获得的长度为L的sh矢量之间的内积得到通过循环右移的时间反转向量x2tr。对于输出样本Yn的1NL时，右循环移位的量为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

22、('Result of circular convolution = ');disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp('Result of IDFT of the DFT products = ');disp(yc)Result of circular convolution =12 28 14 0 16 14Result of IDFT of the DFT products =12 28 14 0 16 14一个圆周卷积的DTF是DTF的逐点产物。Q3.38 运行程序P3

23、.10并验证线性卷积可通过圆周卷积得到。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);ylin = circonv (g1e,g2e);disp('Linear convolution via circular convolution = ');disp(ylin);y = conv(g1, g2);disp('Direct linear convolution = ');disp(y)Linear convolution v

24、ia circular convolution =2 6 10 15 21 15 7 9 5Direct linear convolution =2 6 10 15 21 15 7 9 5使用圆周卷积确实有可能得到线性卷积Q3.40 编写一个MATLAB程序，对两个序列做离散傅里叶变换，已生成他们的线性卷积。用此程序验证Q3.38和Q3.39的结果编写的MATLAB程序：% Program Q3.40g1 = 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 =

25、 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.0000Q3.46 使用程序P3.1在单位圆上求下面的z变换： G(z)=Q3.47 编写一个MATLAB程序，计算并显示零点和极点，计算并显示其因式形式，并产生z的两个多项式之比的形式表示的z变

26、换的极零点图。使用该程序，分析式(3.32)的z变换G(z)。 编写的MATLAB程序：% Program Q3_47clf;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.')

27、;zplane(z,p);运行结果：Zeros: -1.0000 + 1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614iPoles: -8.9576 -0.2718 0.1147 + 0.2627i 0.1147 - 0.2627isos = 1.0000 2.0000 3.0000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822k =0.4000Q3.48 通过习题Q3.47产生的极零点图，求出G(z)的收敛域的数目。清楚地显示所有的收敛域。

28、由极零点图说明离散时间傅里叶变换是否存在。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的存在，它是一个具有双面冲

29、激响应的稳定系统。Q3.50 编写一个MATLAB程序，计算一个有理逆z变换的前L个样本，其中L的值由用户通过命令input提供。用该程序计算并画出式(3.32)中G(z)的逆变换的前50个样本。使用命令stem画出由逆变换产生的序列。编写的MATLAB程序：% Program Q3.50clf;% initializenum = 2 5 9 5 3;den = 5 45 2 1 1;% Query user for parameter LL = input('Enter the number of samples L: ');% find impulse responseg

30、t = impz(num,den,L);%plot the impulse responsestem(t,g);title('First ',num2str(L),' samples of impulse response');xlabel('Time Index n');ylabel('hn');Enter the number of samples L: 50实验名称：数字滤波器的频域分析和实现实验目的和任务：（1）求滤波器的幅度响应和相位响应，观察对称性，判断滤波器类型。（2）用Matlab函数实现系统的级联型和并联型结构，

31、并对滤波器进行结构仿真。实验内容：系统传递函数的级联和并联实现：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('分子系数向量=');den=input('分母系数向量=');z,p,k=tf2zp(num,den);sos=zp2sos(z,p,k);分子系数向量=2 1

32、0 23 34 31 4分母系数向量=1 1 1 1 1 1sos = 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)=分子系数向量 = 3 8 12 7 2 2分母系数向量= 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.25

33、00 1.0000 1.0000 1.0000 1.0000 0.5000 0.5000画出级联实现的框图：Q6.5 使用程序P6.2，式H(z)=所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图：并联 I 形 留数是 -0.4219 + 0.6201i -0.4219 - 0.6201i 2.3437 0.3438 - 2.5079i 0.3438 + 2.5079i极点在 -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i常数 -2 并联II型留数是 -0.30

34、47 - 0.4341i -0.3047 + 0.4341i -1.1719 1.0000 + 0.7758i 1.0000 - 0.7758i极点在 -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i常数 0.1875Q8.1 程序P8.1设计的滤波器类型为 IIR带阻滤波器 ，其指标为 =0.4 0.5 =0.1 0.8 , R=1dB，R=30dB，滤波器的阶数为 4 。为了验证仿真需计算的冲激响应样本为 5个。仿真是正确的吗？Q8.3 生成Q8.1中产生的传输函数的一个级联实现，

35、并编写一个程序来仿真它。其中每个单独的部分用一个直接型实现。验证仿真。format shortnum = input('Numerator coefficients = ');den = input('Denominator coefficients = ');Numfactors = factorize(num);Denfactors = factorize(den);disp('Numerator Factors'),disp(Numfactors)disp('Denominator Factors'),disp(Denfac

36、tors)Numerator coefficients = 0.0571 0 -0.1143 0 0.0571Denominator coefficients = 1.0000 -0.5099 1.2862 -0.3350 0.4479Numerator Factors Columns 1 through 3 1.000000000000000 1.021143238995107 0 1.000000000000000 0.979294541463240 0 1.000000000000000 -1.021143238995070 0 1.000000000000000 -0.97929454

37、1463277 0Denominator Factors Columns 1 through 3 1.000000000000000 -0.597633166080704 0.678485370310279 1.000000000000000 0.087733166080704 0.660146879504817Q8.5 生成Q8.1中传输函数的一个并联型实现，并编写一个程序来仿真它。每一个单独的部分用直接型实现。验证仿真。num = input('Numerator coefficient vector = ');den = input('Denominator co

38、efficient 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);

39、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.0571Denominator coefficient v

40、ector = 1.0000 -0.5099 1.2862 -0.3350 0.4479Parallel Form IResidues are -0.023475623417519 + 0.197736845648730i -0.023475623417519 - 0.197736845648730i -0.011716283258078 - 0.213037612139943i -0.011716283258078 + 0.213037612139943iPoles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352

41、 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311iConstant value 0.127483813351194Parallel Form IIResidues are -0.158795680938490 + 0.041067400713457i -0.158795680938490 - 0.041067400713457i 0.173353325938490 - 0.000160296180734i 0.173353325938490 +

42、 0.000160296180734iPoles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311iConstant value 0.057100000000000b1= -0.046951246835039 -0.289531739579519 0a1= 1.000000000000000 -0.59763316608

43、0704 0.678485370310279Q8.9 程序P8.3设计的滤波器类型为 IIR低通椭圆滤波器 ，指标为=0.25 =0.55，R=0.5dB，R=50dB，阶数为 8，输入的正弦序列的频率为f=0.7kHz。Q8.10 运行P8.3并产生两个图形。哪种输入成分会在滤波器输出出现？为什么输出序列的开始部分不是一种理想的正弦曲线？修改P8.3，以便纸过滤序列X2n。产生的输出序列和预料的是一样的吗？答：产生的输出序列和预料的是不一样的Q8.14 程序P8.4设计的滤波器类型为 FIR低通滤波器，指标为=0 0.3 =0.5 ，阶数为 9 ，为了验证仿真需计算的冲激响应样本为 10个。

44、仿真是正确的吗？实验名称：数字滤波器的设计实验目的和任务：（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

45、定义计算幅度响应 通过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) 根据上述数据

46、和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.2000Q7.5 通过运行P7.1设计巴特沃兹带阻滤波器，写出所产生的传输函数的准确表达式。滤

47、波器的指标是什么？使用MATLAB计算并绘制滤波器未畸变的相位相应及群延迟相应。MATLAB程序为：% Program P7_1Ws = 0.4 0.6; Wp = 0.2 0.8; Rp = 0.4; Rs = 50;N1, Wn1 = buttord(0.2, 0.4, 0.5, 40); num,den = butter(N1,Wn1,'stop'); disp('Numerator Coefficients are ');disp(num);disp('Denominator Coefficients are ');disp(den);%

48、 Compute the gain responseg, w = gain(num,den); % Plot the gain responseplot(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 areColumns 1 through 90.0493 0.0000 0.2465 0.0

49、000 0.4930 0.0000 0.4930 0.0000 0.2465Columns 10 through 110.0 0.0493Denominator Coefficients areColumns 1 through 91.0000 0.0000 -0.0850 0.0000 0.6360 0.0000 -0.0288 0.0000 0.0561Columns 10 through 110.0000 -0.0008因此表达式为：H (z)=指标为=0.2，=0.4，=0.6，=0.8，R=0.4dB，R=50dBQ7.6 修改P7.1后来设计符合习题Q7.1所给指标的切比雪夫1型低

50、通滤波器。写出产生的传输函数的准确表达式。使用MATLAB计算并绘制滤波器未畸变的相位相应及群延迟响应。修改P7.1后的程序（切比雪夫1型低通滤波器）：% Program Q7_6% Design spec as given in Q7.1.FT = 40*103; % sampling freqFp = 4*103; % analog passband edge freqFs = 8*103; % analog stopband edge freqRp = 0.5; % max passband ripple, dBRs = 40; % min stopband attenuation, d

51、B% Convert spec to normalized digital frequenciesomega_p = 2*pi*Fp/FT;Wp = 2*Fp/FT; % omega_p/piomega_s = 2*pi*Fs/FT;Ws = 2*Fs/FT; % omega_s/pi% Estimate the Filter OrderN, Wn = cheb1ord(Wp, Ws, Rp, Rs);% Design the Filternum,den = cheby1(N,Rp,Wn);% Display the transfer functiondisp('Numerator C

52、oefficients are ');disp(num);disp('Denominator Coefficients are ');disp(den);% Compute the gain responseg, w = gain(num,den);% Plot the gain responsefigure(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

53、Chebyshev Lowpass Filter');% Find and plot the phasefigure(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'

54、);% Find and plot the group delayfigure(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_6Ws = 0.4 0.6; Wp = 0.2 0.8; Rp = 0

55、.4; Rs = 50;% Estimate the Filter OrderN1, Wn1 = buttord(Wp, Ws, Rp, Rs);% Design the Filternum,den = butter(N1,Wn1,'stop');% Find the frequency response; find and plot unwrapped phasewp = 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)');t