南京邮电大学DSP实验报告.docx

《南京邮电大学DSP实验报告.docx》由会员分享，可在线阅读，更多相关《南京邮电大学DSP实验报告.docx（34页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、南京邮电高校实 验 报 告试验名称：离散时间信号与系统的时, 频域表示离散傅立叶变换与z变换 数字滤波器的频域分析与实现 数字滤波器的设计 课程名称： 数字信号处理A(双语) 班级学号： 姓 名： 开课时间 ： 2015 /2016 学年 第一学期 试验一试验名称：离散时间信号与系统的时, 频域表示试验目的：熟悉Matlab基本叮嘱与信号处理工具箱，加深理解与驾驭离散时间信号与系统的时, 频域表示及简洁应用。试验任务：在Matlab环境中，依据要求产生序列，对序列进行基本运算；对简洁离散时间系统进行仿真，计算线性时不变（LTI）系统的冲激响应与卷积输出；计算与视察序列的离散时间傅立叶变换（DT

2、FT）幅度谱与相位谱。试验内容：基本序列产生与运算： 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程序：clfn=-10:20;u=zeros(1,10) 1 zeros(1,20);stem(n,u);xlabel(时间序列 n);ylabel(振幅);title(单位样本序列);axis(-10 20 0 1.2);显示的波形如下:Q1.2clf：清除图形；axis：设置坐标轴范围, 可读比例等；title：给图形加标题；xlable：给x轴加标注；

3、ylable：给y轴加标注。Q1.3程序：clfn=-10:20;u=zeros(1,10) 1 zeros(1,20);stem(n+11,u);xlabel(时间序列 n);ylabel(振幅);title(单位样本序列);axis(0 32 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;

4、显示的波形如下:Q1.30 加性噪声dn 是匀整分布在 -0.4与+0.4之间的随机序列Q1.31 不能。因为d是列向量，s是行向量Q1.32 x1是x的延时，x2与x相等，x3超前于xQ1.33 legend用于产生图例说明Q1.30未污染的信号sn 是什么样的形式？加性噪声dn 是什么样的形式？答:未污染的信号sn:是线性增加伴随着实指数缓慢衰减的图像加性噪声dn: 在-0.4与+0.4间匀整分布的自由序列Q1.31运用语句s=s+d能产生被噪声污染的信号吗？若不能，为什么？答:不能，因为- d是一个列向量,而s是一个行向量,须要在添加它们之前调换其中一个向量。Q1.32信号x1, x2,

5、 x3与x之间的关系是什么？答:这三个信号x1,x2,与x3是x扩展的版本,左右边各一个附加的采样。x1是x延迟的版本,一个样本转移到右边并且左边补零。信号x2等于x,左右补0来填充多余的长度。最终,x3是x时间提前的版本,转移一个样本到右边，左边补0。Q1.33legend的作用是什么答:thelegend叮嘱的目的创建图表的说明。在P1_5,信号绘制运用不同的颜色与线类型;说明哪种颜色信息与行类型与每个信号相关联。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(滤波器所需的长度=

6、);滤波器所需的长度=2 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(输入信号);

7、 subplot(2,2,4);plot(n,y);axis(0,100,-2,2);xlabel(时间序列 n);ylabel(振幅);title(输出信号); axis;显示的波形如下:Q2.2程序：n = 0:100;s1 = cos(2*pi*0.05*n); s2 = cos(2*pi*0.47*n); x = s1+s2;M = input(滤波器所需长度 = );num = (-1).0:M-1;y = filter(num,1,x)/M;clf;subplot(2,2,1);plot(n, s1);axis(0, 100, -2, 2);xlabel(时间序号n); ylabe

8、l(振幅);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;显示的波形如下:变更LTI系统对输入的

9、影响是,系统现在是一个高通滤波器。它通过高频输入组件s2来替代低频输入组件s1.Q2.3当M取15时，图像如下Q2.19程序：clf;N = 40;num = 2.2403 2.4908 2.2403;den = 1 -0.4 0.75;y = impz(num,den,N);stem(y);xlabel(时间序号 n); ylabel(振幅);title(冲激响应); grid;显示的波形如下:Q2.21程序：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 = filt

10、er(num,den,x);stem(y);xlabel(时间序号 n); ylabel(振幅);title(冲激响应); grid;显示的波形如下: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(

11、时间序号 n); ylabel(振幅);title(用滤波得到的输出); grid;显示的波形如下:Q3.1答:计算离散时间傅里叶变换的原始序列为：pause叮嘱作用：不加参数,干脆用pause的话,就是程序暂停,直至用户按随意一个按键。假如加参数，例如pause(1)，是程序暂停1秒。Q3.2程序：clf;w = -4*pi:8*pi/511:4*pi;num = 2 1;den = 1 -0.6;h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,real(h);gridtitle(H(ejomega)的实部)xlabel(omega /pi);y

12、label(振幅);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(Amplitude);subplot(2,1,2)plot(w/pi,angle(h);gridtitle( argH(ejomega)的相位谱)xlabel(omega /pi);ylabel(以弧度为单位的相位);显示的波形如

13、下:是w的周期周期是2实部是2为周期是偶对称的；虚部是2为周期是奇对称的；幅度是2为周期是偶对称的；相位是2为周期是奇对称的。Q3.4程序：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( H(ejomega)的实部)xlabel(omega /pi);ylabel(振幅);subplot(2,1,2)plot(w/pi,imag(h);gridtitle( H(ejomega)的虚部)

14、xlabel(omega /pi);ylabel(振幅);pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle( |H(ejomega)|幅度谱)xlabel(omega /pi);ylabel(振幅);subplot(2,1,2)plot(w/pi,angle(h);gridtitle( argH(ejomega)的相位谱)xlabel(omega /pi);ylabel(以弧度为单位的相位);显示的波形如下:试验参考书：S.K.Mitra（著）,孙洪（译）.数字信号处理试验指导书（MATLAB版）. 北京：电子工业出版社，2005试验二试验名称：离散傅

15、立叶变换与z变换试验目的：驾驭离散傅立叶变换（DFT）及逆变换（IDFT）, z变换及逆变换的计算与分析。试验任务：完成DFT与IDFT的计算及常用性质的验证，利用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。试验内容：DFT与IDFT计算： Q3.233.24 （Q3.24可选做）DFT的性质： Q3.263.29，Q3.303.35，Q3.36（Q3.37可选），Q3.38（Q3.39可选），Q3.40z变换分析：Q3.463.48逆z变换：Q3.50试验过程描述：Q3.23程序：clf;N=200; L=256; nn = 0:N-1;kk = 0:L-1;xR = 0.1

16、*(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(实xn);xlabel(时间序号 n);ylabel(振幅);subplot(3,2,2);plot(nn,xI);grid;title(虚xn);xlabel(时间序号 n);ylabel(振幅);subplot(3,2,3);plot(kk,real(XF);grid;title(实Xk);xlabel(频率指数 k);ylabel(振幅);subplot(3,2,4

17、);plot(kk,imag(XF);grid;title(虚Xk);xlabel(频率指数 k);ylabel(振幅);xx = ifft(XF,L);subplot(3,2,5);plot(kk,real(xx);grid;title(IDFTXk实部);xlabel(时间序号 n);ylabel(振幅);subplot(3,2,6);plot(kk,imag(xx);grid;title(IDFTXk虚部);xlabel(时间序号 n);ylabel(振幅);显示的波形如下:Q3.26在函数circshift中，叮嘱rem的作用是什么？答:R=rem(X,Y)，求余数函数，X,Y应当为正

18、数Q3.27说明函数circshift怎样实现圆周移位运算。答:输入序列x是循环左移M位。假如M 0,那么circshift删除左边的元素向量x，并且附加他们到剩下的元素右边来获得循环转移序列。假如假如M 0,然后circshift首先补充的x的长度,最右边的长度(x)- m样品从x中移走并且附加在剩下的M样本右边来得到循环转移序列。 Q3.28在函数circshift中，运算符=的作用是什么？答:假如A与B不相等返回值1 假如A与B相等返回值0Q3.29说明函数circonv怎样实现圆周卷积运算。答:函数circonv操作如下:输入的是两个相等长度为L的两个向量x1与x2.，为了理解circ

19、onv是如何工作的,从x2的周期延拓角度来考虑很有用。让x2p作为x2的无限长的周期延拓。从概念上讲,常规时间反转x2p并且让x2tr 通过x2p的时间反转等于元素1。输出向量y元素1到L是通过x1与一个长度L的通过循环右移一个时间反转序列x2tr得到的序列sh之间的内积来获得的。对于输出样例yn,1nL, 正确的循环移位是n - 1点。Q3.30程序：clf;M = 6;a = 0 1 2 3 4 5 6 7 8 9;b = circshift(a,M); L = length(a)-1;n = 0:L;subplot(2,1,1);stem(n,a);axis(0,L,min(a),max

20、(a);title(原始序列);xlabel(时间序号 n);ylabel(an);subplot(2,1,2);stem(n,b);axis(0,L,min(a),max(a);title(通过循环位移得到的序列 ,num2str(M),样本);xlabel(时间序号 n);ylabel(bn); 确定时移的数量的部分是M假如时移的数量大于序列长度，实际实现的循环时移是rem(M,length(a)点左移，相当于循环移动的M点(不止一次),也相当于通过M点周期延拓的左移。 Q3.31上题程序结果图：序列的长度是10，并且M = 12。这可能被说明为一个12点的循环左移(不止一次),作为一个2

21、点循环左移,或者作为一个2的线性左移，或者序列的12点周期延拓。Q3.32程序：clf;x = 0 2 4 6 8 10 12 14 16; N = length(x)-1; n = 0:N;y = circshift(x,5);XF = fft(x); YF = fft(y); subplot(2,2,1);stem(n,abs(XF);grid;title(原序列的DFT的幅度);xlabel(频率序号 k);ylabel(|Xk|);subplot(2,2,2);stem(n,abs(YF);grid;title(圆周位移后序列的DFT幅度);xlabel(频率序号 k);ylabel(

22、|Yk|);subplot(2,2,3);stem(n,angle(XF);grid;title(原序列的DFT的幅度);xlabel(频率序号 k);ylabel(arg(Xk);subplot(2,2,4);stem(n,angle(YF);grid;title(圆周位移后序列的DFT相位);xlabel(频率序号 k);ylabel(arg(Yk);Q3.33上题程序运行结果图：序列的长度N = 8并且时移是五个样品提前转移到左边。相位是。这是一个重大转变，大大增加了相位谱的斜率。而最初的相位函数只有一个分支切割,在时移信号的相位谱有五个分支切割。Q3.34上述程序修改M=2图像图如下：

23、上述程序修改M=-2结果图如下：Q3.35序列为长度14，上述程序结果图如下：序列长度为16，上述程序结果图如下：Q3.36程序：g1 = 1 2 3 4 5 6; g2 = 1 -2 3 3 -2 1;ycir = cconv(g1,g2);disp(循环卷积图像 = );disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp(DFT变换乘积的IDFT变换的图像 = );disp(yc)结果：循环卷积图像 = Columns 1 through 10 1.0000 0 2.0000 7.0000 10.0000 14.

24、0000 11.0000 28.0000 12.0000 -7.0000 Column 11 6.0000DFT变换乘积的IDFT变换的结果 = 12 28 14 0 16 14Q3.38程序：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 = cconv(g1e,g2e);disp(通过圆周卷积的线性卷积 = );disp(ylin);y = conv(g1, g2);disp(干脆线性卷积 = );disp(y)结果：通过圆周卷积的线性卷积 =

25、Columns 1 through 10 2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.0000 9.0000 5.0000 0.0000 Columns 11 through 17 0.0000 0.0000 0.0000 0 0 0.0000 -0.0000干脆线性卷积 = 2 6 10 15 21 15 7 9 5视察可得: 零填充适当的长度的确可以实现用循环卷积实现线性卷积。Q3.40程序g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros

26、(1,length(g1)-1);G1EF = fft(g1e);G2EF = fft(g2e);ylin = real(ifft(G1EF.*G2EF);disp(通过DFT的线性卷积 = );disp(ylin);结果：通过DFT的线性卷积 = 2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.0000 9.0000 5.0000Q3.46程序：clf;w = 0:pi/51:pi;num = 2 5 9 5 3;den = 5 45 2 1 1;h = freqz(num, den,w);subplot(2,1,1)plot(w/pi,rea

27、l(h);gridtitle(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(振幅);subplot(2,1,2)plot(w/pi,angle(h);gridtitle( argH(ejomega)的相位谱)x

28、label(omega /pi);ylabel(以弧度为单位的相位);图像：Q3.47程序：clf;num = 2 5 9 5 3;den = 5 45 2 1 1;z p k = tf2zpk(num,den);disp(零点:);disp(z);disp(极点:);disp(p);input(按 以接着.);sos k = zp2sos(z,p,k)input(按 以接着.);zplane(z,p);结果：零点: -1.0000 + 1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614i极点: -8.9576 -0.271

29、8 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.0822k = 0.4000按 以接着.Q3.48收敛域数目：R1 : | z | 0.2718 (左边, 不稳定)R 2: 0.2718 | z | 0.2866 (双边, 不稳定)R3 : 0.2866 | z | 8.9576 (右边, 不稳定) 从零极点图可以看出DTFT仅仅从零极点不能得知DTFT是否存在。为了知道,必需指定收敛域。运

30、用上述R3收敛域的序列的DTFT存在。这将是一个稳定的系统且带有一个双边脉冲响应。Q3.50程序：clf;num = 2 5 9 5 3;den = 5 45 2 1 1;L = input(敲入序列长度 L: );g t = impz(num,den,L);stem(t,g);title(前 ,num2str(L), 冲激响应序列);xlabel(时间序号 n);ylabel(hn);图像：试验参考书：S.K.Mitra（著）,孙洪（译）.数字信号处理试验指导书（MATLAB版）. 北京：电子工业出版社，2005试验三试验名称：数字滤波器的频域分析与实现试验目的：驾驭滤波器的传输函数与频率响

31、应的关系，能够从频率响应与零极点模式分析滤波器特性。驾驭滤波器的常用结构。试验任务：求滤波器的幅度响应与相位响应，视察对称性，推断滤波器类型，推断稳定性。验证FIR线性相位滤波器的特点。实现数字滤波器的干脆型, 级联型与并联型结构。试验内容：传输函数与频率响应, 滤波器稳定性：Q4.14.3，Q4.5，Q4.6，Q4.19线性相位滤波器：Q4.19（若群延迟概念没讲，可不求群延迟）数字滤波器结构：Q6.1，Q6.3，Q6.5试验过程描述：Q4.1程序：clear;M = input(键入长度M: );w = 0:2*pi/1023:2*pi;num = (1/M)*ones(1,M);den

32、= 1;h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,abs(h);gridtitle(幅度谱 |H(ejomega)|)xlabel(omega /pi);ylabel(振幅);subplot(2,1,2)plot(w/pi,angle(h);gridtitle(相位谱 argH(ejomega)xlabel(omega /pi);ylabel(以弧度为单位的相位);结果： M=3 M=7： M=10:是低通滤波器Q4.2程序：clf;w = 0:pi/51:pi;num = 0.15 0 -0.15;den = 1 -0.5 0.7;h =

33、freqz(num, den,w);subplot(2,1,1)plot(w/pi,real(h);gridtitle(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(振幅);subplot(2,1,2)plot

34、(w/pi,angle(h);gridtitle( argH(ejomega)的相位谱)xlabel(omega /pi);ylabel(以弧度为单位的相位);图像：带通型滤波器Q4.3程序：clf;w = 0:pi/51:pi;num = 0.15 0 -0.15;den = 0.7 -0.5 1;h = freqz(num, den,w);subplot(2,1,1)plot(w/pi,real(h);gridtitle(H(ejomega)的实部)xlabel(omega /pi);ylabel(振幅);subplot(2,1,2)plot(w/pi,imag(h);gridtitle(

35、 H(ejomega)的虚部)xlabel(omega /pi);ylabel(振幅);pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle( |H(ejomega)|的幅度谱)xlabel(omega /pi);ylabel(振幅);subplot(2,1,2)plot(w/pi,angle(h);gridtitle( argH(ejomega)的相位谱)xlabel(omega /pi);ylabel(以弧度为单位的相位);图像：带通型滤波器与上题相比，本题的滤波器更好Q4.5程序（1）：clf;num = 0.15 0 0.15;den = 1 -0

36、.5 0.7;L = input(敲入序列长度 L: );g t = impz(num,den,L);stem(t,g);title(前 ,num2str(L), 冲激响应序列);xlabel(时间序号 n);ylabel(hn);图像（1）：程序（2）：clf;num = 0.15 0 0.15;den = 0.7 -0.5 1;L = input(敲入序列长度 L: );g t = impz(num,den,L);stem(t,g);title(前 ,num2str(L), 冲激响应序列);xlabel(时间序号 n);ylabel(hn);图像（2）：Q4.6程序（1）：clf;num

37、= 0.15 0 0.15;den = 1 -0.5 0.7;g t = impz(num,den,100);zplane(z,p);图像（1）：程序（2）：clf;num = 0.15 0 0.15;den = 0.7 -0.5 1;g t = impz(num,den,100);zplane(z,p);图像（2）：Q4.19程序:clf;b = 1 -8.5 30.5 -63;num1 = b 81 fliplr(b);num2 = b 81 81 fliplr(b);num3 = b 0 -fliplr(b);num4 = b 81 -81 -fliplr(b);n1 = 0:lengt

38、h(num1)-1;n2 = 0:length(num2)-1;subplot(2,2,1); stem(n1,num1);xlabel(时间序号 n);ylabel(振幅); grid;title(1型有限冲激响应滤波器);subplot(2,2,2); stem(n2,num2);xlabel(时间序号 n);ylabel(振幅); grid;title(2型有限冲激响应滤波器);subplot(2,2,3); stem(n1,num3);xlabel(时间序号 n);ylabel(振幅); grid;title(3型有限冲激响应滤波器);subplot(2,2,4); stem(n2,n

39、um4);xlabel(时间序号 n);ylabel(振幅); grid;title(4型有限冲激响应滤波器);pausesubplot(2,2,1); zplane(num1,1);title(1型有限冲激响应滤波器);subplot(2,2,2); zplane(num2,1);title(2型有限冲激响应滤波器);subplot(2,2,3); zplane(num3,1);title(3型有限冲激响应滤波器);subplot(2,2,4); zplane(num4,1);title(4型有限冲激响应滤波器);disp(1型有限冲激响应滤波器的零点是);disp(roots(num1);

40、disp(2型有限冲激响应滤波器的零点是);disp(roots(num2);disp(3型有限冲激响应滤波器的零点是);disp(roots(num3);disp(4型有限冲激响应滤波器的零点是);disp(roots(num4);结果：1型有限冲激响应滤波器的零点是 2.9744 2.0888 0.9790 + 1.4110i 0.9790 - 1.4110i 0.3319 + 0.4784i 0.3319 - 0.4784i 0.4787 0.3362 2型有限冲激响应滤波器的零点是 3.7585 + 1.5147i 3.7585 - 1.5147i 0.6733 + 2.6623i 0.6733 - 2.6623i -1.0000 0.0893 + 0.3530i 0.