《南京邮电大学DSP试验报告.doc》由会员分享，可在线阅读，更多相关《南京邮电大学DSP试验报告.doc（75页珍藏版）》请在三一文库上搜索。
1、南京邮电大学实验报告实验名称： 离散时间信号与系统的时、频域表示离散傅立叶变换和 z变换数字滤波器的频域分析和实现数字滤波器的设计课程名称：数字信号处理A（双语）班级学号:姓 名:开课时间：2015 /2016 学年 第一学期实验一实验名称：离散时间信号与系统的时、频域表示实验目的：熟悉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程序：clfn=-10:20;u=zeros（1,10） 1 zeros（1,20）;stem（ n,u）;xlabel（' 时间序列 n'）;ylabel（'振幅'）;title（' 单位样本序列'）;axis（-10 20 0 1.2）;单位样本序列幅 振时间序列nQ1.2elf :清除图形；axis：设置坐标轴范
3、围、可读比例等;title :给图形加标题；xlable :给x轴加标注；ylable :给y轴加标注。Q1.3程序：elfn=-10:20;u=zeros(1,10) 1 zeros(1,20);stem (n+11,u);xlabel('时间序列n');ylabel('振幅');title('单位样本序列');axis(0 32 0 1.2);显示的波形如下10.8幅060.4 -单位样本序列0510152025时间序列n30Q1.23程序：n=0:50;>> f=0.08;>> phase=pi/2;>>
4、 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;>>显示的波形如下正弦序列幅 振时间序列nQ1.30加性噪声dn是均匀分布在 -0.4和+0.4之间的随机序
5、列Q1.31不能。因为d是列向量，s是行向量Q1.32x1是x的延时，x2和x相等，x3超前于xQ1.33lege nd用于产生图例说明Q1.30未污染的信号sn是什么样的形式？加性噪声dn是什么样的形式?答:未污染的信号sn:是线性增加伴随着实指数缓慢衰减的图像加性噪声dn:在-0.4和+0.4间均匀分布的自由序列Q1.31使用语句s=s+d能产生被噪声污染的信号吗？若不能，为什么？答:不能，因为-d是一个列向量，而s是一个行向量，需要在添加它们之前调换其中一个向量。Q1.32信号x1、x2、x3与x之间的关系是什么？答:这三个信号x1,x2,和x3是x扩展的版本，左右边各一个附加的采样。x
6、1是x延迟的版本,x3一个样本转移到右边并且左边补零。信号x2等于x,左右补0来填充多余的长度。最后是x时间提前的版本，转移一个样本到右边，左边补 0。Q1.33lege nd的作用是什么答:thelegend命令的目的创建图表的说明。在P1_5,信号绘制使用不同的颜色和线类型说明哪种颜色信息和行类型与每个信号相关联。Q2.1程序：elf;>> n=0:100;>> s仁 cos(2*pi*0.05* n);>> s2=cos(2*pi*0.47* n);>> x=s1+s2;>> M=input('滤波器所需的长度='
7、;);滤波器所需的长度=2>> num=on es(1,M);>> y=filter( nu m,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);&
8、gt;> 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('振幅'); tit
9、le('输出信号');>> axis;显示的波形如下2-2信号#110-1050100时间序列n输入信号信号#2时间序列n050100时间序列n幅 振Q2.2程序：n = 0:100;s1 = cos(2*pi*0.05* n);s2 = cos(2*pi*0.47* n);x = s1+s2;M = in put('滤波器所需长度=');num = (-1).a0:M-1;y = filter( num,1,x)/M;clf;subplot(2,2,1);plot (n, s1);axis(0, 100, -2, 2);振幅');xlabe
10、l('时间序号 n'); ylabel('ti tle('信号 #1');subplot(2,2,2);plot (n, s2);axis(0, 100, -2, 2);振幅');xlabel(' 时间序号 n'); ylabel('ti tle('信号 #2');subplot(2,2,3);plot( n, x);axis(0, 100, -2, 2);xlabel(' 时间序号 n'); ylabel('振幅');title(' 输入信号');subpl
11、ot(2,2,4);Plot( n, y)；axis(0, 100, -2, 2);振幅');xlabel(' 时间序号 n'); ylabel(' title('输岀信号');axis;显示的波形如下：幅 振210-1-221幅 0振-1-2信号#2210-1-2050100时间序号n 输入信号050100时间序号n 输出信号21幅°振-1-2050100时间序号n050100时间序号ns2来替代改变LTI系统对输入的影响是，系统现在是一个高通滤波器。它通过高频输入组件 低频输入组件s1.Q2.3信号#121:i!0-i-2100幅
12、振50时间序号n幅 振信号#22时间序号n 输出信号050100时间序号n当M取15时，图像如下10-1-2Q2.19程序：elf;N = 40;num = 2.2403 2.4908 2.2403;den = 1 -0.4 0.75;y = impz（ num,de n, N）;stem（y）;xlabel（'时间序号 n' ）; ylabel（ '振幅）; title（'冲激响应'）;grid;显示的波形如下冲激响应幅 振时间序号nQ2.21程序:elf；N = 40;num = 0.9 -0.45 0.35 0.002;den = 1.0 0.71
13、 -0.46 -0.62;x = 1 zeros(1,N-1);y = filter( num,de n,x);stem(y);振幅');xlabel(' 时间序号 n'); ylabel(' title('冲激响应');grid;显示的波形如下幅 振冲激响应时间序号nQ2.28程序：elf;h = 3 2 1 -2 1 0 -4 0 3;x = 1 -2 3 -4 3 2 1;y = con v(h,x);n = 0:14;subplot(2,1,1);stem( n, y);xlabel(' 时间序号 n'); ylabel(
14、'振幅');ti tle('用卷积得到的输岀');grid;x1 = x zeros(1,8);y1 = filter(h,1,x1);subplot(2,1,2);stem( n, y1);xlabel(' 时间序号 n'); ylabel('振幅');ti tle('用滤波得到的输岀');grid;显示的波形如下20用卷积得到的输出100-10-2014200-10-201012146 8时间序号niiA10126 8时间序号n用滤波得到的输出丄A10Q3.1答:计算离散时间傅里叶变换的原始序列为:H(ejw)
15、2 z1 -0.6zJpause命令作用：不加参数，直接用pause的话,就是程序暂停，直至用户按任意一个按键。如 果加参数，例如pause(1),是程序暂停1秒。Q3.2 程序: elf;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);grid title('H(eyomega)的实部')xlabel('omega 八pi');ylabel(' 振幅');subplot(2,1,2)的虚
16、部')plot(w/pi,imag(h);grid ti tle(' H(eAjomega) xlabel('omega 八pi'); ylabel(' 振幅'); pause subplot(2,1,1) plot(w/pi,abs(h);grid title(' |H(eAjomega)|的幅度谱')xlabel('omega 八pi');ylabel('Amplitude'); subplot(2,1,2) plot(w/pi,a ngle(h);grid title(' argH(e
17、Ajomega)的相位谱')xlabel('omega 八pi');ylabel('以弧度为单位的相位');显示的波形如下幅振H(e J')的实部- / 二H(e j )的虚部幅|H(e )|的幅度谱-10 12- / :argH(eh-.)的相位谱34位相的位单为度弧以210-1-2-4-3-2-10 12- / 二34是W的周期周期是2 n实部是2 n为周期是偶对称的; 虚部是2n为周期是奇对称的; 幅度是2 n为周期是偶对称的;相位是2n为周期是奇对称的。Q3.4程序：elf；w = -4*pi:8*pi/511:4*pi;num = 1
18、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(' H(eAjomega) xlabel('omega 八pi');ylabel(' 振幅'); subplot(2,1,2) plot(w/pi,imag(h);grid title(' H(eAjomega) xlabel('omega 八pi');ylabel(' 振幅');pause subplot(2,1,1)
19、plot(w/pi,abs(h);grid title(' |H(eAjomega)|的实部')的虚部')幅度谱')xlabel('omega 八pi'); ylabel(' 振幅'); subplot(2,1,2) plot(w/pi,a ngle(h);grid title(' argH(eAjomega)的相位谱')xlabel('omega 八pi');ylabel('以弧度为单位的相位');显示的波形如下100H（e）的实部100500-500- / 二H(ej )的虚部5
20、01a111| A«j;A Cl-;f|川1 Vl V 1/ 1 1 r1JTTVr y'1 ; ；!11 1 I1IJ,1/IILV0-50-4-3-2A(I卩l 1 : A A rA1 :. / /A A A 11 J V v :.1 v ! 1"1| V 八1'八'Jfifu-4-3-20- /-100ir/1.111111f1111X jJ11111iIJ-|H(ej )|幅度谱幅50-3-2-1 0 1- / 二argH(e j )的相位谱位相的位单为度弧以1 11 11 1 1Il i 1LI | 1 1 '- 1N11!iHi&
21、#39;.l1iI''!. I11 |'1JniW)il;|'：,1,put 1H4 2 0-2-3-2实验参考书：SKMitra （著），孙洪（译）数字信号处理实验指导书（ MATLAB版）.北京：电子工 业出版社，2005实验二实验名称：离散傅立叶变换和z变换实验目的：掌握离散傅立叶变换（DFT）及逆变换（IDFT ）、z变换及逆变 换的计算和分析。实验任务：完成DFT和IDFT的计算及常用性质的验证，利用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。实验内容：DFT 和 IDFT 计算：Q3.233.24 （Q3.24可选做）DFT 的性质：
22、Q3.263.29, Q3.303.35, Q3.36 （Q3.37 可选），Q3.38 （Q3.39 可选），Q3.40z变换分析：Q3.463.48逆z变换：Q3.50实验过程描述：Q3.23程序：elf;N=200;L=256;nn = 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('实 x n'); xlabel('时间序号n&
23、#39;); ylabel('振幅')； subplot(3,2,2); plot (nn, xI);grid; title('虚 x n'); xlabel('时间序号n'); ylabel('振幅')； subplot(3,2,3); plot(kk,real(XF);grid; title('实 Xk'); xlabel('频率指数k'); ylabel('振幅')； subplot(3,2,4); plot(kk,imag(XF);grid; title('虚 Xk&
24、#39;); 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('振幅')；显示
25、的波形如下实xn实Xk幅 振虚xn虚Xk幅 振10000-1000300幅 振100 200频率指数kIDFTXk实部500-500rL 一1r'¥I-00300100-10300100 200时间序号n时间序号n幅 振100 200频率指数kIDFTXk虚部Q3.26在函数circshift中，命令rem的作用是什么？ 答:R=rem(X,Y)，求余数函数，X,Y应该为正数Q3.27解释函数circshift怎样实现圆周移位运算。答:输入序列x是循环左移 M位。如果M > 0,那么circshift删除左边的元素向量X，并且附加他们到剩下的元素右边来获得循环转移序列。
26、如果如果M < 0,然后circshift首先补充的x的长度，最右边的长度(x)- m样品从x中移走并且附加在剩下的 M样本右边来得到循环转移 序列。Q3.28在函数circshift中，运算符=的作用是什么？答:如果A和B不相等返回值1如果A和B相等返回值0Q3.29解释函数circonv怎样实现圆周卷积运算。答:函数circonv操作如下：输入的是两个相等长度为L的两个向量x1和x2.，为了理解circonv是如何工作的，从x2的周期延拓角度来考虑很有用。让 x2p作为x2的无限长的周期延拓。从概念上讲，常规时间反转 x2p并且让x2tr通过x2p的时间反转等于元素1。输出向量y元素
27、1到L是通过x1和一个长度L的通过循环右移一个时间反转序列x2tr得到的序列sh之间的内积来获得的。对于输出样例 yn,1 w nw L、正确的循环移位是n - 1点。Q3.30程序：elf;M = 6;a = 0 1 2 3 4 5 6 7 8 9;b = circshift(a,M);L = len gth(a)-1;n = 0:L;subplot(2,1,1);stem( n,a);axis(O,L,mi n(a),max(a);title('原始序列)；xlabel('时间序号n');ylabel('a n');subplot(2,1,2);st
28、em( n,b);axis(O,L,mi n(a),max(a);title('通过循环位移得到的序列,num2str(M),'样本');xlabel('时间序号n');ylabel('b n');决定时移的数量的部分是 M如果时移的数量大于序列长度，实际实现的循环时移是rem(M,length(a)点左移，相当于循环移动的M点(不止一次),也相当于通过 M点周期延拓的左移。Q3.31上题程序结果图:原始序列通过循环位移得到的序列6样本序列的长度是10 ,并且M = 12。这可能被解释为一个12点的循环左移(不止一次),作为一个2点循环左
29、移，或者作为一个2的线性左移，或者序列的12点周期延拓。Q3.32程序：elf;x = 0 2 4 6 8 10 12 14 16;N = len gth(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('圆周位移
30、后序列的DFT幅度');xlabel(' 频率序号k');ylabel('|Yk|');subplot(2,2,3);stem (n,an gle(XF);grid;title('原序列的DFT的幅度');xlabel(' 频率序号k');ylabel('arg(Xk)');subplot(2,2,4);stem (n,an gle(YF);grid;title('圆周位移后序列的DFT相位');xlabel(' 频率序号k');ylabel('arg(丫k)'
31、;);Q3.33上题程序运行结果图:WNkn)8060X 4020原序列的DFT的幅度G)(080C60丫 4020圆周位移后序列的DFT幅度 T <1002468频率序号k原序列的DFT的幅度02468频率序号k圆周位移后序列的DFT相位的长度N = 8并且时移是个样品提前转移到左边。相位是_k 5叫 =ejk10 7/8jk5 血 4 =e 。这是一个重大转变，大大增加了相位谱的斜率。而最初的相位函数只有一个分支切割，在时移信号的相位谱有五个分支切割。Q3.34上述程序修改M=2图像图如下:8060X 4020原序列的DFT的幅度G)(080C60Y 4020圆周位移后序列的DFT幅
32、度 T <1002468频率序号k原序列的DFT的幅度02468频率序号k圆周位移后序列的DFT相位频率序号k上述程序修改M=-2结果图如下:原序列的DFT的幅度圆周位移后序列的DFT幅度o o o O8 6 4 2806020、1K246频率序号k原序列的DFT的幅度00246频率序号k4 2 0-2圆周位移后序列的DFT相位2468频率序号k-40246频率序号kQ3.35序列为长度14,上述程序结果图如下:2001009.!tiI99I?r；'：15015200500510频率序号k原序列的DFT的幅度A序列长度为16,上述程序结果图如下:X r(j)11 11?r&quo
33、t;15010050051015频率序号k200kX10051015频率序号kmm300200100原序列的DFT的幅度频率序号k1015频率序号kQ3.36程序：g1 = 1 2 3 4 5 6; g2 = 1 -2 3 3 -2 1;ycir = ccon v(g1,g2);disp('循环卷积图像=');disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp('DFT变换乘积的IDFT变换的图像=');disp(yc)结果: 循环卷积图像 =Colu mns 1 through 10
34、1.000012.0000-7.00000 2.00007.000010.000014.000011.000028.0000Colu mn 116.0000DFT变换乘积的IDFT变换的结果=12281401614Q3.38程序：g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,le ngth(g2)-1);g2e = g2 zeros(1,le ngth(gl)-l);yli n = ccon v(g1e,g2e);disp('通过圆周卷积的线性卷积=');disp(yli n);y = con v(g1, g2);disp(
35、9;直接线性卷积 =');disp(y)结果：通过圆周卷积的线性卷积=Colu mns 1 through 102.00006.000010.000015.000021.000015.00007.00005.00000.00009.0000Colu mns 11 through 170.00000.00000.0000000.0000 -0.0000直接线性卷积=2 610152115795观察可得：零填充适当的长度确实可以实现用循环卷积实现线性卷积。Q3.40程序g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,le ngth(g2)-1)
36、;g2e = g2 zeros(1,le ngth(g1)-1);G1EF = fft(g1e);G2EF = fft(g2e);ylin = real(ifft(G1EF.*G2EF);disp('通过 DFT 的线性卷积=');disp(ylin);结果：通过DFT的线性卷积 =9.00002.00006.000010.000015.000021.000015.00007.00005.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
37、(2,1,1)plot(w/pi,real(h);gridtitle('H(eyomega)的实部')xlabel('omega 八pi');ylabel(' 振幅');subplot(2,1,2)plot(w/pi,imag(h);gridtitle(' H(eAjomega)的虚部')xlabel('omega 八pi');ylabel(' 振幅');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle(' |H(eAjomega)|的幅度谱'
38、;)xlabel('omega 八pi');ylabel(' 振幅');subplot(2,1,2)plot(w/pi,a ngle(h);gridtitle(' argH(eAjomega)的相位谱')xlabel('omega 八pi');ylabel('以弧度为单位的相位');图像:00.60.40.2-0.20 0.1H(ej )的实部0.20.30.40.50.60.7- / 二H(ej )的虚部0.80.9-0.1-0.2-0.31-0.400.10.20.30.40.50.60.70.80.9- /
39、二0.80.6幅04 振0.40.2|H(ej )|的幅度谱000.10.20.30.40.50.60.70.80.9- / 二argH(e j )的相位谱0位相的位单为度弧以-1-2-3-400.10.20.30.40.50.60.70.80.9- / 二Q3.47程序:Cif;num = 2 5 9 5 3;den = 5 45 2 1 1;z p k=tf2zpk (n um,de n);disp('零点:');disp(z);disp('极点:');disp(p);in put('按 <return>以继续sos k=:zp2sos(
40、z,p,k)in put('按 <return>以继续zpla ne(z,p);结果：零点：-1.0000 + 1.4142i-1.0000 - 1.4142i-0.2500 + 0.6614i-0.2500 - 0.6614i极点：-8.9576-0.27180.1147 + 0.2627i0.1147 - 0.2627i按<return>以继续sos1.00002.00003.00001.00009.22932.43441.00000.50000.50001.0000-0.22930.0822k =0.4000按return以继续Q3.48收敛域数目：Ri
41、: | z | < 0.2718左边,不稳定)R 2： 0.2718 < | z | < 0.2866 双边,不稳定)R3： 0.2866 < | z | < 8.9576 双边,稳定)R4： | z | > 8.9576 (右边,不稳定)从零极点图可以看出DTFT 仅仅从零极点不能得知DTFT是否存在。为了知道必须指定收敛域。使用上述R3收敛域的序列的 DTFT存在。这将是一个稳定的系统且带有一个双边脉冲响应。Q3.50程序：elf；num = 2 5 9 5 3; den = 5 45 2 1 1;L = in put('敲入序列长度 L:');g t = impz( num,de n,L);stem(t,g);title('前',num2str(L),'冲激响应序列');xlabel('时间序号n');ylabel('h n');