您的当前位置：首页南邮DSP实验报告

南邮DSP实验报告

来源：意榕旅游网

2014-2015第二学期DSP上机实验报告

南京邮电大学

实验报告

实验名称：离散时间信号与系统的时、频域表示

离散傅立叶变换和z变换

数字滤波器的频域分析和实现

数字滤波器的设计

课程名称：数字信号处理A(双语)

班级学号： B12020

姓名：

开课时间： 2014 /2015 学年第二学期

2014-2015第二学期DSP上机实验报告

实验一

实验名称：离散时间信号与系统的时、频域表示实验目的：

熟悉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 程序： 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]); 显示的波形如下:

2014-2015第二学期DSP上机实验报告

单位样本序列10.8振幅0.60.40.20-10-505时间序列 n101520

Q1.2

clf：清除图形；

axis：设置坐标轴范围、可读比例等； title：给图形加标题； xlable：给x轴加标注； ylable：给y轴加标注。

Q1.3 程序： clf

n=-10:20;

u=[zeros(1,10) 1 zeros(1,20)]; stem(n+11,u);

xlabel('时间序列 n');ylabel('振幅'); title('单位样本序列'); axis([0 32 0 1.2]);

显示的波形如下:

2014-2015第二学期DSP上机实验报告

单位样本序列10.8振幅0.60.40.2005101520时间序列 n2530

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; >>

显示的波形如下:

2014-2015第二学期DSP上机实验报告

正弦序列321振幅0-1-2-3051015202530时间序列n304550

Q1.30

加性噪声d[n] 是均匀分布在 -0.4和+0.4之间的随机序列

Q1.31

不能。因为d是列向量，s是行向量

Q1.32

x1是x的延时，x2和x相等，x3超前于x

Q1.33

legend用于产生图例说明

Q1.30

未污染的信号s[n] 是什么样的形式？加性噪声d[n] 是什么样的形式？答:未污染的信号s[n]:是线性增加伴随着实指数缓慢衰减的图像加性噪声d[n]: 在-0.4和+0.4间均匀分布的自由序列

Q1.31

使用语句s=s+d能产生被噪声污染的信号吗？若不能，为什么？

答:不能，因为- d是一个列向量,而s是一个行向量,需要在添加它们之前调换其中一个向量。

Q1.32

信号x1、x2、x3与x之间的关系是什么？

答:这三个信号x1,x2,和x3是x扩展的版本,左右边各一个附加的采样。x1是x延迟的版本,

2014-2015第二学期DSP上机实验报告

一个样本转移到右边并且左边补零。信号x2等于x,左右补0来填充多余的长度。最后,x3是x时间提前的版本,转移一个样本到右边，左边补0。

Q1.33

legend的作用是什么

答: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('滤波器所需的长度='); 滤波器所需的长度=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('输入信号'); >> subplot(2,2,4); plot(n,y);

axis([0,100,-2,2]);

xlabel('时间序列 n');ylabel('振幅'); title('输出信号'); >> axis;

显示的波形如下:

2014-2015第二学期DSP上机实验报告

信号#12121信号#2振幅0-1-2050时间序列 n输入信号100振幅0-1-2050时间序列 n输出信号1002121振幅0-1-2050时间序列 n100振幅0-1-2050时间序列 n100

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'); 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]);

2014-2015第二学期DSP上机实验报告

xlabel('时间序号 n'); ylabel('振幅'); title('输入信号'); subplot(2,2,4); plot(n, y);

axis([0, 100, -2, 2]);

xlabel('时间序号 n'); ylabel('振幅'); title('输出信号'); axis;

显示的波形如下:

信号 #121振幅信号 #221振幅0-1-2050时间序号n输入信号1000-1-2050 时间序号n输出信号10021振幅21振幅0-1-2050时间序号 n1000-1-2050时间序号 n100

改变LTI系统对输入的影响是,系统现在是一个高通滤波器。它通过高频输入组件s2来替代低频输入组件s1.

Q2.3

2014-2015第二学期DSP上机实验报告

当M取15时，图像如下

信号 #121振幅信号 #221振幅0-1-2050时间序号n输入信号1000-1-2050 时间序号n输出信号10021振幅21振幅0-1-2050时间序号 n1000-1-2050时间序号 n100

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;

2014-2015第二学期DSP上机实验报告

显示的波形如下:

冲激响应4321振幅0-1-2-30510152025时间序号 n3030

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 = filter(num,den,x); stem(y);

xlabel('时间序号 n'); ylabel('振幅'); title('冲激响应'); grid;

2014-2015第二学期DSP上机实验报告

显示的波形如下:

冲激响应21.510.5振幅0-0.5-1-1.50510152025时间序号 n3030

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;

2014-2015第二学期DSP上机实验报告

显示的波形如下:

用卷积得到的输出2010振幅0-10-2002468时间序号 n用滤波得到的输出1012142010振幅0-10-2002468时间序号 n101214

Q3.1

2z1答:计算离散时间傅里叶变换的原始序列为：H(e) 110.6zjw

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));grid title('H(e^{j\\omega})的实部') xlabel('\\omega /\\pi'); ylabel('振幅'); subplot(2,1,2)

2014-2015第二学期DSP上机实验报告

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('Amplitude'); subplot(2,1,2)

plot(w/pi,angle(h));grid

title(' arg[H(e^{j\\omega})]的相位谱') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

显示的波形如下:

H(ej)的实部86振幅420-4-3-2-10 / H(ej)的虚部421234振幅0-2-4-4-3-2-10 /1234

2014-2015第二学期DSP上机实验报告

|H(ej)|的幅度谱8Amplitude20-4-3-2-101234 / arg[H(ej)]的相位谱210-1-2-4-3-2-101234以弧度为单位的相位 /

是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));grid title(' H(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})|幅度谱')

2014-2015第二学期DSP上机实验报告

xlabel('\\omega /\\pi'); ylabel('振幅'); subplot(2,1,2)

plot(w/pi,angle(h));grid

title(' arg[H(e^{j\\omega})]的相位谱') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

显示的波形如下:

H(ej)的实部10050振幅0-50-4-3-2-10 / H(ej)的虚部123410050振幅0-50-100-4-3-2-10 /1234

2014-2015第二学期DSP上机实验报告

|H(ej)|幅度谱100振幅500-4-3-2-10 /1234 arg[H(ej)]的相位谱以弧度为单位的相位420-2-4-4-3-2-10 /1234

实验参考书：

S.K.Mitra（著）,孙洪（译）.数字信号处理实验指导书（MATLAB版）. 北京：电子工业出版社，2005

2014-2015第二学期DSP上机实验报告

实验二

实验名称：离散傅立叶变换和z变换实验目的：

掌握离散傅立叶变换（DFT）及逆变换（IDFT）、z变换及逆变换的计算和分析。实验任务：

完成DFT和IDFT的计算及常用性质的验证，利用DFT实现线性卷积，实现z变换的零极点分析，求有理逆z变换。实验内容：

DFT和IDFT计算： Q3.23～3.24 （Q3.24可选做）

DFT的性质： Q3.26~3.29，Q3.30~3.35，Q3.36（Q3.37可选），

Q3.38（Q3.39可选），Q3.40

z变换分析：Q3.46～3.48 逆z变换：Q3.50 实验过程描述：

Q3.23

程序：

clf;

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;

2014-2015第二学期DSP上机实验报告

XF = fft(x,L);

subplot(3,2,1);grid; plot(nn,xR);grid; title('实\\{x[n]\\}'); xlabel('时间序号 n'); 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('实\\{X[k]\\}'); xlabel('频率指数 k'); ylabel('振幅'); subplot(3,2,4);

plot(kk,imag(XF));grid; title('虚\\{X[k]\\}'); xlabel('频率指数 k'); ylabel('振幅'); xx = ifft(XF,L); subplot(3,2,5);

plot(kk,real(xx));grid; title('IDFT\\{X[k]\\}实部'); xlabel('时间序号 n'); ylabel('振幅'); subplot(3,2,6);

plot(kk,imag(xx));grid; title('IDFT\\{X[k]\\}虚部'); xlabel('时间序号 n'); ylabel('振幅');

2014-2015第二学期DSP上机实验报告

显示的波形如下:

实{x[n]}101虚{x[n]}振幅50050100150时间序号 n实{X[k]}200振幅0-1050100150时间序号 n虚{X[k]}2001000500振幅0-10000100200频率指数 kIDFT{X[k]}实部300振幅0-5000100200频率指数 kIDFT{X[k]}虚部300101振幅0-100100200时间序号 n300振幅0-10100200时间序号 n300

Q3.26

在函数circshift中，命令rem的作用是什么？答:R=rem(X,Y)，求余数函数，X,Y应该为正数 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相等返回值0

Q3.29

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

答:函数circonv操作如下:输入的是两个相等长度为L的两个向量x1和x2.，为了理解circonv是如何工作的,从x2的周期延拓角度来考虑很有用。让x2p作为x2的无限长的周期延拓。

2014-2015第二学期DSP上机实验报告

从概念上讲,常规时间反转x2p并且让x2tr 通过x2p的时间反转等于元素1。输出向量y元素1到L是通过x1和一个长度L的通过循环右移一个时间反转序列x2tr得到的序列sh之间的内积来获得的。对于输出样例y[n],1≤n≤L、正确的循环移位是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(a)]); title('原始序列'); xlabel('时间序号 n'); ylabel('a[n]'); subplot(2,1,2);

stem(n,b);axis([0,L,min(a),max(a)]);

title(['通过循环位移得到的序列 ',num2str(M),'样本']); xlabel('时间序号 n'); ylabel('b[n]');

决定时移的数量的部分是M

如果时移的数量大于序列长度，实际实现的循环时移是rem(M,length(a))点左移，相当于循环移动的M点(不止一次),也相当于通过M点周期延拓的左移。

Q3.31

上题程序结果图：

2014-2015第二学期DSP上机实验报告

原始序列86a[n]420012456时间序号 n通过循环位移得到的序列 6样本3786b[n]420012345时间序号 n67

序列的长度是10，并且M = 12。这可能被解释为一个12点的循环左移(不止一次),作为一个2点循环左移,或者作为一个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('|X[k]|'); subplot(2,2,2); stem(n,abs(YF));grid;

title('圆周位移后序列的DFT幅度'); xlabel('频率序号 k'); ylabel('|Y[k]|'); subplot(2,2,3);

stem(n,angle(XF));grid; title('原序列的DFT的幅度'); xlabel('频率序号 k'); ylabel('arg(X[k])'); subplot(2,2,4);

2014-2015第二学期DSP上机实验报告

stem(n,angle(YF));grid;

title('圆周位移后序列的DFT相位'); xlabel('频率序号 k'); ylabel('arg(Y[k])');

Q3.33

上题程序运行结果图：

原序列的DFT的幅度80608060圆周位移后序列的DFT幅度|X[k]|40200046频率序号 k原序列的DFT的幅度28|Y[k]|402000468频率序号 k圆周位移后序列的DFT相位24242arg(X[k])0-2-40246频率序号 k8arg(Y[k])0-2-40246频率序号 k8

序列的长度N = 8并且时移是五个样品提前转移到左边。相位是

kn0k5WNWNejk10/8ejk5/4。

这是一个重大转变，大大增加了相位谱的斜率。而最初的相位函数只有一个分支切割,在时移信号的相位谱有五个分支切割。

2014-2015第二学期DSP上机实验报告

Q3.34

上述程序修改M=2图像图如下：

上述程序修改M=-2结果图如下：

2014-2015第二学期DSP上机实验报告

Q3.35

序列为长度14，上述程序结果图如下：

2014-2015第二学期DSP上机实验报告

原序列的DFT的幅度200150200150圆周位移后序列的DFT幅度|X[k]|100|Y[k]|100500051015频率序号 k原序列的DFT的幅度42)]k[X(0rga-2-4051015频率序号 k

序列长度为16，上述程序结果图如下：

500051015频率序号 k圆周位移后序列的DFT相位420-2-4051015频率序号 karg(Y[k])

2014-2015第二学期DSP上机实验报告

原序列的DFT的幅度300200300200圆周位移后序列的DFT幅度|X[k]|1000|Y[k]|10000510频率序号 k原序列的DFT的幅度15051015频率序号 k圆周位移后序列的DFT相位4242arg(X[k])0-2-40510频率序号 k15arg(Y[k])0-2-40510频率序号 k15

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.0000 11.0000 28.0000 12.0000 -7.0000

Column 11

6.0000

DFT变换乘积的IDFT变换的结果 =

12 28 14 0 16 14 Q3.38

2014-2015第二学期DSP上机实验报告

程序：

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)

结果：

通过圆周卷积的线性卷积 = Columns 1 through 10

2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.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(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 5.0000

Q3.46 程序：

clf;

w = 0:pi/51:pi;

9.0000 9.0000 2014-2015第二学期DSP上机实验报告

num = [2 5 9 5 3];den = [5 45 2 1 1]; h = freqz(num, den,w); subplot(2,1,1)

plot(w/pi,real(h));grid title('H(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(' arg[H(e^{j\\omega})]的相位谱') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

图像：

2014-2015第二学期DSP上机实验报告

H(ej)的实部0.60.4振幅0.20-0.200.10.20.30.40.5 /0.60.70.80.91 H(ej)的虚部0-0.1振幅-0.2-0.3-0.400.10.20.30.40.5 /0.60.70.80.91

|H(ej)|的幅度谱0.80.6振幅0.40.2000.10.20.30.40.5 /0.60.70.80.91 arg[H(ej)]的相位谱以弧度为单位的相位0-1-2-3-400.10.20.30.40.5 /0.60.70.80.91

Q3.47

2014-2015第二学期DSP上机实验报告

程序：

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.2718 0.1147 + 0.2627i 0.1147 - 0.2627i

按以继续...

sos =

1.0000 2.0000 3.0000 1.0000 9.2293 1.0000 0.5000 0.5000 1.0000 -0.2293 k =

0.4000

按以继续...

2.4344 0.0822 2014-2015第二学期DSP上机实验报告

432Imaginary Part10-1-2-3-4-9-8-7-6-5-4-3Real Part-2-101

Q3.48

收敛域数目：

R1 : | 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 程序：

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');

2014-2015第二学期DSP上机实验报告

ylabel('h[n]');

图像：

20-2-4-6x 1045前 50 冲激响应序列h[n]-8-10-12-14-16051015202530时间序号 n304550

实验参考书：

S.K.Mitra（著）,孙洪（译）.数字信号处理实验指导书（MATLAB版）. 北京：电子工业出版社，2005

2014-2015第二学期DSP上机实验报告

实验三

实验名称：数字滤波器的频域分析和实现实验目的：

掌握滤波器的传输函数和频率响应的关系，能够从频率响应和零极点模式分析滤波器特性。掌握滤波器的常用结构。实验任务：

求滤波器的幅度响应和相位响应，观察对称性，判断滤波器类型，判断稳定性。验证FIR线性相位滤波器的特点。实现数字滤波器的直接型、级联型和并联型结构。实验内容：

传输函数和频率响应、滤波器稳定性：Q4.1～4.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 = [1];

h = freqz(num, den, w); subplot(2,1,1)

plot(w/pi,abs(h));grid

2014-2015第二学期DSP上机实验报告

title('幅度谱 |H(e^{j\\omega})|') xlabel('\\omega /\\pi'); ylabel('振幅'); subplot(2,1,2)

plot(w/pi,angle(h));grid

title('相位谱 arg[H(e^{j\\omega})]') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

结果：

M=3

幅度谱 |H(ej)|1振幅0.5000.20.40.60.81 /1.21.41.61.82相位谱 arg[H(ej)]以弧度为单位的相位420-2-400.20.40.60.81 /1.21.41.61.82

M=7：

2014-2015第二学期DSP上机实验报告

幅度谱 |H(ej)|1振幅0.5000.20.40.60.81 /1.21.41.61.82相位谱 arg[H(ej)]以弧度为单位的相位420-2-400.20.40.60.81 /1.21.41.61.82

M=10:

幅度谱 |H(ej)|1振幅0.5000.20.40.60.81 /1.21.41.61.82相位谱 arg[H(ej)]以弧度为单位的相位420-2-400.20.40.60.81 /1.21.41.61.82

是低通滤波器

2014-2015第二学期DSP上机实验报告

Q4.2 程序：

clf;

w = 0:pi/51:pi;

num = [0.15 0 -0.15];den = [1 -0.5 0.7]; h = freqz(num, den,w); subplot(2,1,1)

plot(w/pi,real(h));grid title('H(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(' arg[H(e^{j\\omega})]的相位谱') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

2014-2015第二学期DSP上机实验报告

图像：

H(ej)的实部1振幅0.5000.10.20.30.40.5 /0.60.70.80.91 H(ej)的虚部0.5振幅0-0.500.10.20.30.40.5 /0.60.70.80.91

|H(ej)|的幅度谱1振幅0.5000.10.20.30.40.50.60.70.80.91 / arg[H(ej)]的相位谱210-1-200.10.20.30.40.5 /0.60.70.80.91以弧度为单位的相位

带通型滤波器

Q4.3 程序：

clf;

w = 0:pi/51:pi;

num = [0.15 0 -0.15];den = [0.7 -0.5 1];

2014-2015第二学期DSP上机实验报告

h = freqz(num, den,w); subplot(2,1,1)

plot(w/pi,real(h));grid title('H(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(' arg[H(e^{j\\omega})]的相位谱') xlabel('\\omega /\\pi'); ylabel('以弧度为单位的相位');

图像：

H(ej)的实部0.50振幅-0.5-100.10.20.30.40.5 /0.60.70.80.91 H(ej)的虚部0.5振幅0-0.500.10.20.30.40.50.60.70.80.91 /

2014-2015第二学期DSP上机实验报告

|H(ej)|的幅度谱1振幅0.5000.10.20.30.40.50.60.70.80.91 / arg[H(ej)]的相位谱420-2-400.10.20.30.40.5 /0.60.70.80.91以弧度为单位的相位

带通型滤波器

与上题相比，本题的滤波器更好

Q4.5

程序（1）：

clf;

num = [0.15 0 0.15]; den = [1 -0.5 0.7];

L = input('敲入序列长度 L: '); [g t] = impz(num,den,L); stem(t,g);

title(['前 ',num2str(L),' 冲激响应序列']); xlabel('时间序号 n'); ylabel('h[n]');

2014-2015第二学期DSP上机实验报告

图像（1）：

前 100 冲激响应序列0.150.10.05h[n]0-0.05-0.10102030405060时间序号 n708090100

程序（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('h[n]');

2014-2015第二学期DSP上机实验报告

图像（2）：

3210x 106前 100 冲激响应序列h[n]-1-2-3-4-50102030405060时间序号 n708090100

Q4.6

程序（1）：

clf;

num = [0.15 0 0.15]; den = [1 -0.5 0.7]; [g t] = impz(num,den,100); zplane(z,p);

2014-2015第二学期DSP上机实验报告

图像（1）：

10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.51

程序（2）：

clf;

num = [0.15 0 0.15]; den = [0.7 -0.5 1]; [g t] = impz(num,den,100); zplane(z,p);

2014-2015第二学期DSP上机实验报告

图像（2）：

10.80.60.4Imaginary Part0.20-0.2-0.4-0.6-0.8-1-1-0.50Real Part0.51

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:length(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,num4);

xlabel('时间序号 n');ylabel('振幅'); grid;

2014-2015第二学期DSP上机实验报告

title('4型有限冲激响应滤波器'); pause

subplot(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));

disp('2型有限冲激响应滤波器的零点是'); disp(roots(num2));

disp('3型有限冲激响应滤波器的零点是'); disp(roots(num3));

disp('4型有限冲激响应滤波器的零点是'); disp(roots(num4));

结果：

1型有限冲激响应滤波器10050100502型有限冲激响应滤波器振幅0-50-1000468时间序号 n3型有限冲激响应滤波器2振幅0-50-1000510时间序号 n4型有限冲激响应滤波器1005010050振幅0-50-1000246时间序号 n8振幅0-50-10005时间序号 n10

2014-2015第二学期DSP上机实验报告

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.03 + 0.3530i 0.03 - 0.3530i 0.22 + 0.0922i 0.22 - 0.0922i

3型有限冲激响应滤波器的零点是 4.7627 1.6279 + 3.0565i 1.6279 - 3.0565i -1.0000 1.0000 0.1357 + 0.29i 0.1357 - 0.29i 0.2100

4型有限冲激响应滤波器的零点是 3.4139 1.61 + 1.5813i 1.61 - 1.5813i -0.0733 + 0.9973i -0.0733 - 0.9973i 1.0000 0.3159 + 0.3020i 0.3159 - 0.3020i 0.2929

2014-2015第二学期DSP上机实验报告

1型有限冲激响应滤波器10-1-1123Real Part3型有限冲激响应滤波器10-108210-1-22型有限冲激响应滤波器Imaginary PartImaginary Part9-224Real Part4型有限冲激响应滤波器0Imaginary Part0-2-28Imaginary Part2902Real Part46-1012Real Part3

Q6.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)

2014-2015第二学期DSP上机实验报告

结果： 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 得到：

H1(z)2(13z12z2)(1z12z2)(1z10.5z2)

Q6.3 程序：

num = [3 8 12 7 2 -2]; den = [16 24 24 14 5 1];

[b,a] = eqtflength(num,den); % make lengths equal [z,p,k] = tf2zp(b,a); sos = zp2sos(z,p,k)

结果： 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

2014-2015第二学期DSP上机实验报告

Q6.5

程序：

num = [3 8 12 7 2 -2]; den = [16 24 24 14 5 1];

[r1,p1,k1] = residuez(num,den); [r2,p2,k2] = residue(num,den); disp('并联I型')

disp('留数是');disp(r1); disp('极点在');disp(p1); disp('常数');disp(k1); disp('并联II型')

disp('留数是');disp(r2); disp('极点在');disp(p2); disp('常数');disp(k2);

2014-2015第二学期DSP上机实验报告

结果：并联I型留数是

-0.4219 + 0.6201i -0.4219 - 0.6201i 2.3437 0.3437 - 2.5079i 0.3437 + 2.5079i

极点在

-0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i 常数 -2

并联II型留数是

-0.3047 - 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.1875

框图为：

2014-2015第二学期DSP上机实验报告

实验参考书：

S.K.Mitra（著）,孙洪（译）.数字信号处理实验指导书（MATLAB版）. 北京：电子工业出版社，2005

2014-2015第二学期DSP上机实验报告

实验四

实验名称：数字滤波器设计

实验目的：掌握数字滤波器的常用设计方法，在Matlab软件环境下设计满足给定指标的滤波器，综合运用滤波器的不同设计手段和滤波器的应用知识，设计出符合一定应用要求的滤波器。

实验任务：以模拟低通滤波器为原型设计IIR数字滤波器，包括根据指标选择阶数、确定滤波器的传输函数。用窗口法设计FIR数字滤波器，包括观察吉布斯现象。选定一个滤波问题，设计数字滤波器，呈现滤波效果。实验内容：

（1）IIR数字滤波器设计：Q7.1（阶数估计），Q7.5～7.6（滤波器设计）

（2）FIR数字滤波器设计：Q7.13和/或Q7.14（阶数估计），Q7.20（滤波器设计）（3）数字滤波器应用

含噪序列为x[n]=s[n]+w[n]，其中s[n]为幅度等于5的直流信号,w[n]是方差为１的高斯白噪声。要求对500点序列进行去噪运算。请自行选择滤波器类型，确定滤波器截止频率和阶数并设计滤波器。画出滤波器的频率响应，画出输入和输出信号的波形和频谱，计算降

2014-2015第二学期DSP上机实验报告

噪前后的信噪比增益。实验过程描述：

Q7.1

最后计算得 FT40KHz FP4KHz FS8KHz

RP0.5dB RS40dB

Q7.5 程序：

Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50; [N1, Wn1] = buttord(Wp, Ws, Rp, Rs); [num,den] = butter(N1,Wn1,'stop'); disp('分子系数是 ');disp(num); disp('分母系数是 ');disp(den); [g, w] = gain(num,den); plot(w/pi,g);grid axis([0 1 -60 5]);

xlabel('\\omega /\\pi'); ylabel('增益， dB'); title('巴特沃兹带组滤波器的增益响应');

结果：

分子系数是

Columns 1 through 10

0.0493 0.0000 0.2465 0.0000 0.4930 0.0000 0.2465 0.0000

Column 11

0.0493

分母系数是

Columns 1 through 10

1.0000 0.0000 -0.0850 0.0000 0.6360 0.0000 0.0561 0.0000

Column 11

-0.0008

0.4930 -0.0288 0.0000 0.0000 2014-2015第二学期DSP上机实验报告

Q7.6

程序：

FT = 40*10^3; Fp = 4*10^3; Fs = 8*10^3; Rp = 0.5; Rs = 40;

omega_p = 2*pi*Fp/FT; Wp = 2*Fp/FT;

omega_s = 2*pi*Fs/FT; Ws = 2*Fs/FT;

[N, Wn] = cheb1ord(Wp, Ws, Rp, Rs); [num,den] = cheby1(N,Rp,Wn); disp('分子系数是 ');disp(num); disp('分母系数是 ');disp(den); [g, w] = gain(num,den); figure(1);

plot(w/pi,g);grid; axis([0 1 -60 5]);

xlabel('\\omega /\\pi'); ylabel('增益， dB'); title(' 1 型Chebyshev低通滤波器增益响应'); % 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('非约束相位 (rad)'); title(' 1 型Chebyshev 低通滤波器非约束相位响应'); % Find and plot the group delay figure(3);

GR = grpdelay(num,den,w2); plot(w2/pi,GR);grid;

xlabel('\\omega /\\pi'); ylabel('群延迟 (sec)'); title('1型 Chebyshev 低通滤波器群延迟');

2014-2015第二学期DSP上机实验报告

结果：

分子系数是

0.0004 0.0020 0.0040 0.0040 0.0020 0.0004

分母系数是

1.0000 -3.8269 6.2742 -5.44 2.4915 -0.4797

2014-2015第二学期DSP上机实验报告

Q7.13

经过计算 WP2KHz WS2KHz P0.005 S0.005

FT10KHz N=46

Q7.20 程序：

clear; Fp = 2*10^3; Fs = 2.5*10^3; FT = 10*10^3; Rp = 0.005; Rs = 0.005;

N = kaiord(Fp,Fs,Rp,Rs,FT) Wp = 2*Fp/FT; Ws = 2*Fs/FT;

Wn = Wp + (Ws - Wp)/2; h = fir1(N,Wn);

disp('分子系数是 ');disp(h); [g, w] = gain(h,[1]); figure(1);

plot(w/pi,g);grid;

xlabel('\\omega /\\pi'); ylabel('增益 dB'); title('增益响应'); w2 = 0:pi/511:pi; Hz = freqz(h,[1],w2); MagH = abs(Hz);

T1 = 1.005*ones(1,length(w2)); T2 = 0.995*ones(1,length(w2)); T3 = 0.005*ones(1,length(w2)); figure(4);

plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid; figure(2);

Phase = angle(Hz); plot(w2/pi,Phase);grid;

xlabel('\\omega /\\pi'); ylabel('Phase (rad)'); title('频率响应'); figure(3);

UPhase = unwrap(Phase); plot(w2/pi,UPhase);grid;

xlabel('\\omega /\\pi'); ylabel('非约束相位 (rad)'); title('非约束相位响应');

2014-2015第二学期DSP上机实验报告

结果：

0.0010 -0.0004 -0.0015 0.0000 0.0024 0.0010 -0.0038 -0.0032 0.0049 0.0071 -0.0050 -0.0128 0.0026 0.0202 0.0038 -0.0284 -0.0166 0.0366 0.0404 -0.0436 -0.0909 0.0483 0.3129 0.4498 0.3129 0.0483 -0.0909 -0.0436 0.0404 0.0366 -0.0166 -0.0284 0.0038 0.0202 0.0026 -0.0128 -0.0050 0.0071 0.0049 -0.0032 -0.0038 0.0010 0.0024 0.0000 -0.0015 -0.0004 0.0010

2014-2015第二学期DSP上机实验报告

阶数为66

实验参考书：

S.K.Mitra（著）,孙洪（译）.数字信号处理实验指导书（MATLAB版）. 北京：电子工业出版社，2005

2014-2015第二学期DSP上机实验报告

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文