DSP实验报告封面及要求.doc

《DSP实验报告封面及要求.doc》由会员分享，可在线阅读，更多相关《DSP实验报告封面及要求.doc（56页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、南京邮电大学实 验 报 告实验名称：离散时间信号与系统的时、频域表示离散傅立叶变换和z变换 数字滤波器的频域分析和实现数字滤波器的设计课程名称 数字信号处理A(双语) 班级学号_B_姓 名_吴兵_开课时间 2013/2014学年，第2学期 实验名称：离散时间信号与系统的时、频域表示实验目的和任务：熟悉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(Time index n);ylabel(Amplitude);title(Unit Sample Sequence);axis(-10 20 0 1.2);Q1.2 clf 命令的作用是- 清除图形窗口上的图形

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

4、。该序列的周期是多少 n = 0:50; f = 0.08; phase = 90; A = 2.5; arg = 2*pi*f*n - phase; x = A*cos(arg); clf; % Clear old graph stem(n,x); % Plot the generated sequence axis(0 50 -3 3); grid; title(Sinusoidal Sequence); xlabel(Time index n); ylabel(Amplitude); axis; Q1.29 运行程序P1.5的图形Q1.30 未污染的信号sn 是什么样的形式？加性噪声dn

5、 是什么样的形式？ 未污染的信号加性噪声dn 是均匀分布在 -0.4和+0.4之间的随机序列Q1.31 使用语句s=s+d能产生被噪声污染的信号吗？若不能，为什么？不能。因为d是列向量，s是行向量Q1.32 信号x1、x2、x3与x之间的关系是什么？x1是x延时一个单位，x2和x相等，x3超前于x一个单位Q1.33 legend的作用是什么legend用于产生图例说明Q2.1 对于M = 2 和输入 xn = s1n+s2n，程序P2.1的输出为：. 输入 xn 被该离散时间系统抑制的分量为- 高频分量Q2.2 程序P2.1 中 LTI system 被修改为 yn = 0.5(xnxn1)后

6、， 输入 xn = s1n+s2n 导致的输出为:对于输入的影响是-该系统现在是一个高通滤波器。它通过高频率的输入分量，抑制低频输入分Q2.3程序 P2_1 对于不同M（M=5,10）取值和不同正弦分量(任取2个)取值的运行结果如下： M=5 f1=0.01 f=0.5 M=10 f1=0.05 f2=0.9Q2.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,de

7、n,x);stem(y);xlabel(时间序号n); ylabel(振幅);title(冲击响应); grid;程序产生的40个样本如下所示: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计算

8、离散时间傅里叶变换的原始序列为- pause 命令的作用为- 程序暂停,直至用户按任意一个按键,程序才继续运行Q3.2 程序 P3_1 运行结果为:DTFT 是关于 w的周期函数么？-DTFT 是关于 w的周期函数周期是 - 四个图形的对称性为： 实部是偶对称，虚部是奇对称，幅度谱是偶对称，相位谱是奇对称。Q3.4 修改程序 P3_1 重做Q3.2的程序如下:clf;w = -4*pi:8*pi/511:4*pi;num = 1 3 5 7 9 11 13 1 17;den = 1;h = freqz(num, den, w);% Plot the DTFTsubplot(2,1,1)plot

9、(w/pi,real(h);gridtitle(Real part of H(ejomega)xlabel(omega /pi);ylabel(Amplitude);subplot(2,1,2)plot(w/pi,imag(h);gridtitle(Imaginary part of H(ejomega)xlabel(omega /pi);ylabel(Amplitude);pause subplot(2,1,1)plot(w/pi,abs(h);gridtitle(Magnitude Spectrum |H(ejomega)|)xlabel(omega /pi);ylabel(Amplitu

10、de);subplot(2,1,2)plot(w/pi,angle(h);gridtitle(Phase Spectrum argH(ejomega)xlabel(omega /pi);ylabel(Phase in radians);修改程序后的运行结果为: DTFT 是关于 w的周期函数么？-DTFT 是关于 w的周期函数周期是 - 相位谱中跳变的原因是 - 对相位进行了归一化 实验名称：离散傅立叶变换和z变换实验目的和任务：掌握离散傅立叶变换（DFT）及逆变换（IDFT）、z变换及逆变换的计算和分析。利用Matlab语言，完成DFT和IDFT的计算及常用性质的验证，用DFT实现线性卷积，

11、实现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点离散傅里叶逆变换想xn。对不同长度N和不同的离散傅里叶变换长度L，运行程序。讨论你的结果。程序：%N=200,L=256%长度为N的L点DFTclf;N=200; % length of signalL=256; % length of DFTLON =

12、 0:N-1;LOL = 0:L-1;x = 0.1*(1:100) zeros(1,N-100); XF = fft(x,L);subplot(3,2,1);grid;plot(LON,x);grid;title(xn);xlabel(Time index n);ylabel(Amplitude); subplot(3,2,3);plot(LOL,real(XF);grid;title(ReXk);xlabel(Frequency index k);ylabel(Amplitude); subplot(3,2,4);plot(LOL,imag(XF);grid;title(ImXk);xla

13、bel(Frequency index k);ylabel(Amplitude); % IDFTxx = ifft(XF,L);subplot(3,2,5);plot(LOL,real(xx);grid;title(Real part of IDFTXk);xlabel(Time index n);ylabel(Amplitude);subplot(3,2,6);plot(LOL,imag(xx);grid;title(Imag part of IDFTXk);xlabel(Time index n);ylabel(Amplitude);运行结果：Q3.24 写一个MATLAB程序，用一个N点

14、复数离散傅里叶计算两个长度为N的实数序列的N点离散傅里叶变换，并将结果同直接使用两个N点离散傅里叶变换得到的结果进行比较。程序：g=1 2 1.5 1 4 ;h=2.1 2 1 1.5 3;x=1+i*2.1 2+i*2 1.5+i*1 1+i*1.5 4+i*3;G=fft(g);H=fft(h);X=fft(x);for n=1:5;XK(n)=conj(X(mod(-(n-1),5)+1);endG1=0.5*(X+XK);H1=-0.5*i*(X-XK);运行结果：G = 9.5000 + 0.0000i 0.8316 + 1.6082i -3.0816 + 1.6511i -3.08

15、16 - 1.6511i 0.8316 - 1.6082i G1G1 = 9.5000 + 0.0000i 0.8316 + 1.6082i -3.0816 + 1.6511i -3.0816 - 1.6511i 0.8316 - 1.6082i HH = 9.6000 + 0.0000i 1.6225 + 1.2449i -1.1725 + 0.1123i -1.1725 - 0.1123i 1.6225 - 1.2449i H1H1 = 9.6000 + 0.0000i 1.6225 + 1.2449i -1.1725 + 0.1123i -1.1725 - 0.1123i 1.6225

16、- 1.2449i显然，两者的计算结果完全一致。Q3.26 在函数circshift中，命令rem的作用是什么？答：rem(x,y)是用y对x求余数函数。Q3.27 解释函数circshift怎样实现圆周移位运算。答：在输入序列x由M的位置开始被循环移位。如果M 0，则circshift删除从矢量x最左边开始的M个元素和它们附加在右侧的剩余元素，以获得循环移位序列。如果如果M0，则circshift首先通过x的长度来弥补M，即序列x最右边的长度的M样品从x中删除和所附在其余的M个样本的右侧，以获得循环移位序列。Q3.28 在函数circshift中，运算符=的作用是什么？ 答：=是不等于的意思

17、。Q3.29 解释函数circonv怎样实现圆周卷积运算。 答：输入是两个长度都为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 =

18、Circonv(g1,g2);disp(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 14Q3.38 运行程序P3.10并验证线性卷积可通过圆

19、周卷积得到。程序：% Program P3_10% Linear Convolution via Circular Convolutiong1 = 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 = );

20、disp(y)运行结果：Linear convolution via 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的结果程序：g1 = 2 1 3 5 4;g2 = 1 2 3 0 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros(1,length(g1)

21、-1);G1=fft(g1e);G2=fft(g2e);yIDFT=real(ifft(G1.*G2);disp(Linear convolution via IDFT = );disp(yIDFT);y = conv(g1, g2);disp(Direct linear convolution = );disp(y)运行结果：Linear convolution via IDFT = 2.0000 5.0000 11.0000 14.0000 25.0000 24.0000 15.0000 5.0000 4.0000Direct linear convolution = 2 5 11 14

22、25 24 15 5 4再次验证Q3.46 使用程序P3.1在单位圆上求下面的z变换： G(z)=Q3.47 编写一个MATLAB程序，计算并显示零点和极点，计算并显示其因式形式，并产生z的两个多项式之比的形式表示的z变换的极零点图。使用该程序，分析式(3.32)的z变换G(z)。程序：clf;num = 2 5 9 5 3;den = 5 45 2 1 1;zplane(num,den);z p k = tf2zpk(num,den);disp(Zeros:)disp(z)disp(Poles:)disp(p)sos = zp2sos(z,p,k)运行结果：Zeros: -1.0000 +

23、1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614iPoles: -8.9576 + 0.0000i -0.2718 + 0.0000i 0.1147 + 0.2627i 0.1147 - 0.2627isos = 0.4000 0.8000 1.2000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822Q3.48 通过习题Q3.47产生的极零点图，求出G(z)的收敛域的数目。清楚地显示所有的收敛域。由极零点图说明离散时间傅里叶变换是否存在。R1

24、: | z | 0.2718 (left-sided, not stable)R2 : 0.2718 | z | 0.2866 (two-sided, not stable)R 3: 0.2866 | z | 8.9576 (right-sided, not stable) 不能从极零点图肯定的说DTFT是否存在，因为其收敛域一定要指定。当收敛域在上述R 3内所获得的序列却是证明了DTFT的存在，它是一个具有双面冲激响应的稳定系统。Q3.50 编写一个MATLAB程序，计算一个有理逆z变换的前L个样本，其中L的值由用户通过命令input提供。用该程序计算并画出式(3.32)中G(z)的逆变换的

25、前50个样本。使用命令stem画出由逆变换产生的序列。程序：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(hn);实验名称：数字滤波器的频域分析和实现实验目的：掌握滤波器的传输函数和频率响应的关系，能够从频率响应和零极点模式分析滤波器特性。掌握滤波器的常用结

26、构。实验任务：求滤波器的幅度响应和相位响应，观察对称性，判断滤波器类型，判断稳定性。验证FIR线性相位滤波器的特点。实现数字滤波器的直接型、级联型和并联型结构。实验内容：传输函数和频率响应、滤波器稳定性：Q4.14.3，Q4.5，Q4.6，Q4.19线性相位滤波器：Q4.19数字滤波器结构：Q6.1，Q6.3，Q6.5数字滤波器仿真：Q8.1，Q 8.3，Q 8.5，Q 8.9，Q 8.10，Q 8.14实验过程与结果分析：Q4.1 修改程序P2.1,取三个不同的M值，当时计算并画出式（2.13）所示滑动平均滤波器的幅度和相位谱。证明由幅度和相位谱表现出的对称型。它表示了哪种类型的滤波器？修改

27、后的程序：clf;w = 0:8*pi/511:2*pi;M=input(滤波器所需长度=);num = ones(1,M);y = freqz(num, 1, w)/M;subplot(2,1,1)plot(w/pi,abs(y);gridtitle(Magnitude Spectrum |Y(ejomega)|)xlabel(omega /pi);ylabel(Amplitude);subplot(2,1,2)plot(w/pi,angle(y);gridtitle(Phase Spectrum argY(ejomega)xlabel(omega /pi);ylabel(Phase in

28、radians);运行结果： M=2 M=4 M=10 它表示了II型滤波器。Q4.2 使用修改后的程序P3.1,计算并画出当时传输函数 的因果线性时不变离散时间系统的频率响应。它表示哪种类型的滤波器？ 修改后的程序：clf;w = 0:8*pi/511:pi;M=input(滤波器所需长度=);num =0.15 0 -0.15;den=1 -0.5 0.7;y = freqz(num, den, w)/M;subplot(2,1,1)plot(w/pi,abs(y);gridtitle(Magnitude Spectrum |Y(ejomega)|)xlabel(omega /pi);yl

29、abel(Amplitude);subplot(2,1,2)plot(w/pi,angle(y);gridtitle(Phase Spectrum argY(ejomega)xlabel(omega /pi);ylabel(Phase in radians); 运行结果： 它表示带通滤波器。Q4.3 对下面的传输函数重做Q4.2: 修改后的程序：clf;w = 0:8*pi/511:pi;M=input(滤波器所需长度=);num =0.15 0 -0.15;den=0.7 -0.5 1;y = freqz(num, den, w)/M;subplot(2,1,1)plot(w/pi,abs(

30、y);gridtitle(Magnitude Spectrum |Y(ejomega)|)xlabel(omega /pi);ylabel(Amplitude);subplot(2,1,2)plot(w/pi,angle(y);gridtitle(Phase Spectrum argY(ejomega)xlabel(omega /pi);ylabel(Phase in radians);运行结果：区别在于相位谱，我会选择前者来滤波，因为前者的相位谱相对平缓，后者在通带内有较大的相位突变。Q4.5 使用习题Q3.50中编写的程序，分别计算并画出式（4.36）和式（4.37）确定的两个滤波器的冲击

31、响应的前100个样本。讨论你的结果。修改后的程序：式4.36clf;num = 0.15 0 -1;den = 0.7 -0.5 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(hn);运行结果：式4.37修改后的程序：clf;num = 0.15 0 -1;den = 0.7 -0.5 1;L = input(Enter th

32、e 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(hn);运行结果： 两者反转。Q4.6 使用zplane分别生成式（4.36）和式（4.37）确定的两个滤波器的零点图。讨论你的结果。程序：（1）num = 0.15 0 -1;den = 1 -0.5 0.7;zplane(num,den);(2) num = 0.15 0 -1;den = 0.7 -0.5 1;zp

33、lane(num,den);运行结果：（1） （2）两者的零点一样，极点前者在圆内，后者在圆外。Q4.19 运行程序P4.3, 生成每一类线性相位有限冲击响应滤波器的冲击响应。每一个有限冲击响应滤波器的长度是多少？验证冲击序列的对称性。接着验证这些滤波器的零点位置。使用MATLAB计算并绘出这些滤波器的相位响应，验证他们的线性相位特性。这些滤波器的群延迟是多少？运行结果：（1）长度分别是9,10,9,10；（2）Zeros of Type 1 FIR Filter are 2.9744 + 0.0000i 2.0888 + 0.0000i 0.9790 + 1.4110i 0.9790 - 1

34、.4110i 0.3319 + 0.4784i 0.3319 - 0.4784i 0.4787 + 0.0000i 0.3362 + 0.0000iZeros of Type 2 FIR Filter are 3.7585 + 1.5147i 3.7585 - 1.5147i 0.6733 + 2.6623i 0.6733 - 2.6623i -1.0000 + 0.0000i 0.0893 + 0.3530i 0.0893 - 0.3530i 0.2289 + 0.0922i 0.2289 - 0.0922iZeros of Type 3 FIR Filter are 4.7627 + 0.

35、0000i 1.6279 + 3.0565i 1.6279 - 3.0565i -1.0000 + 0.0000i 1.0000 + 0.0000i 0.1357 + 0.2549i 0.1357 - 0.2549i 0.2100 + 0.0000iZeros of Type 4 FIR Filter are 3.4139 + 0.0000i 1.6541 + 1.5813i 1.6541 - 1.5813i -0.0733 + 0.9973i -0.0733 - 0.9973i 1.0000 + 0.0000i 0.3159 + 0.3020i 0.3159 - 0.3020i 0.2929

36、 + 0.0000i（3）程序：pausew = -4*pi:8*pi/511:4*pi;h1=freqz(num1,1,w);h2=freqz(num2,1,w);h3=freqz(num3,1,w);h4=freqz(num4,1,w);subplot(2,2,1); plot(w/pi,angle(h1);gridtitle(Phase Spectrum)xlabel(omega /pi);ylabel(Phase in radians); grid;title(Type 1 FIR Filter);subplot(2,2,2); plot(w/pi,angle(h2);gridtitl

37、e(Phase Spectrum)xlabel(omega /pi);ylabel(Phase in radians); grid;title(Type 2 FIR Filter);subplot(2,2,3); plot(w/pi,angle(h3);gridtitle(Phase Spectrum)xlabel(omega /pi);ylabel(Phase in radians); grid;title(Type 3 FIR Filter);subplot(2,2,4); plot(w/pi,angle(h4);gridtitle(Phase Spectrum)xlabel(omega

38、/pi);ylabel(Phase in radians); grid;title(Type 4 FIR Filter);运行结果：（4）程序：pausesubplot(2,2,1);grpdelay(h1);title(Type 1 FIR Filter);subplot(2,2,2);grpdelay(h1);title(Type 2 FIR Filter);subplot(2,2,3);grpdelay(h1);title(Type 3 FIR Filter);subplot(2,2,4);grpdelay(h1);title(Type 4 FIR Filter);Q6.1 使用程序P6

39、.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 0H(z)不是一个线性相位传输函数。Q6.3 使用程序P6.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)=所示因果无限冲激响应传输函数的两种不同并联形式实现。画出两种实现的框图：并联 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.4330