数字信号处理书上实验1.2.3.4

实验一 熟悉Matlab环境

一、实验目的

1. 熟悉MATLAB的主要操作命令。 2. 学会简单的矩阵输入和数据读写。 3. 掌握简单的绘图命令。

4. 用MATLAB编程并学会创建函数。 5. 观察离散系统的频率响应。 二、实验内容

认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。

上机实验内容：

（1）数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。

clear all; a=[1 2 3 4]; b=[3 4 5 6]; c=a+b; d=a-b; e=a.*b; f=a./b; g=a.^b; n=1:4;

subplot(4,2,1);stem(n,a);

xlabel('n');xlim([0 5]);ylabel('A'); subplot(4,2,2);stem(n,b);

xlabel('n');xlim([0 5]);ylabel('B'); subplot(4,2,3);stem(n,c);

xlabel('n');xlim([0 5]);ylabel('C'); subplot(4,2,4);stem(n,d); xlabel('n');xlim([0 5]);ylabel('D'); subplot(4,2,5);stem(n,e);

xlabel('n');xlim([0 5]);ylabel('E'); subplot(4,2,6);stem(n,f); xlabel('n');xlim([0 5]);ylabel('F'); subplot(4,2,7);stem(n,g); xlabel('n');xlim([0 5]);ylabel('G');

（2）用MATLAB实现下列序列：

a) x(n)=0.8n 0≤n≤15

b) x(n)=e(0.2+3j)n 0≤n≤15

c) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15

d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。

e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。

clear all; N=0:15;

% a) x(n)=0.8n 0≤n≤15 xa=0.8.^N;

- 1 -

figure;subplot(2,1,1);stem(N,xa); xlabel('n');xlim([0 16]);ylabel('xa'); % b) x(n)=e(0.2+3j)n 0≤n≤15 xb=exp((0.2+3*j)*N); subplot(2,1,2);stem(N,xb);

xlabel('n');xlim([0 16]);ylabel('xb');figure;

% c) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15 xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);

subplot(3,1,1);stem(N,xc);xlabel('n');xlim([0 16]);ylabel('xc');

% d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。 k=0:3;m=0; for i=1:4

for j=1:16 m=m+1;

n(m)=N(j)+16*k(i);

x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi); end

end

subplot(3,1,2);stem(n,x16);xlabel('n');ylabel('x16');

% e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。 for j=1:10 x10(j)=x16(j); end

for i=1:3

for m=1:10 x10(i*10+m)=x10(m); end end n=1:40;

subplot(3,1,3);stem(n,x10); xlabel('n');ylabel('x10');

- 2 -

（3）x(n)=[1,-1,3,5],产生并绘出下列序列的样本: a) x1(n)=2x(n+2)-x(n-1)-2x(n)

5b) x2(n)?clear all n=1:4;

?nx(n?k)

k?1T=4;

x=[1 -1 3 5]; x(5:8)=x(1:4);

subplot(2,1,1);stem(1:8,x);grid;

for i=1:4 if i-1<0

x1(i)=2*x(i+2)-x(i-1)-2*x(i); else

x1(i)=2*x(i+2)-x(i-1+T)-2*x(i); end end

x1(5:8)=x1(1:4);

subplot(2,1,2);stem(1:8,x1);grid;

（4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注： a) x(t)=sin(2πt) 0≤t≤10s

b) x(t)=cos(100πt)sin(πt) 0≤t≤4s

ta=0:0.05:10;

xa=sin(2*pi*ta);

subplot(2,1,1);plot(ta,xa); xlabel('t');ylabel('幅度'); tb=0:0.01:4;

- 3 -

xb=cos(100*pi*tb).*sin(pi*tb); subplot(2,1,2);plot(tb,xb); xlabel('t');ylabel('幅度');

（5）编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

终点为n2。

n0=5;ns=1;nf=10;%ns为起点;nf为终点；在=n=n0处生成单位阶跃序列 n=[ns:nf];

x=[(n-n0)>=0]; stem(n,x);

（6）给一定因果系统H(z)与相频响应。

clear all; b=[1,sqrt(2),1]; a=[1,-0.67,0.9]; [h,w]=freqz(b,a);

am=20*log10(abs(h)); subplot(2,1,1);plot(w,am);

?(1?2z-1?1)/(1-0.67z-1?0.9z-2)求出并绘制H(z)的幅频响应

- 4 -

ph=angle(h);

subplot(2,1,2);plot(w,ph);

（7）计算序列{8 -2 -1 2 3}和序列{2 3 -1 -3}的离散卷积，并作图表示卷积结果。

clear all;

a=[8 -2 -1 2 3]; b=[2 3 -1 -3];

c=conv(a,b); %计算卷积 M=length(c)-1; n=0:1:M; stem(n,c);

xlabel('n');ylabel('幅度');

（8）求以下差分方程所描述系统的单位脉冲响应h(n),0≤n≤50 y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1)

clear all; N=50; a=[1 -2];

b=[1 0.1 -0.06]; x=[1 zeros(1,N-1)]; k=0:1:N-1;

y=filter(a,b,x); stem(k,y);

xlabel('n');ylabel('幅度 ');

- 5 -

实验三 IIR数字滤波器的设计

(1)fc?0.3kHz，??0.8dB，fr?0.2kHz，At?20dB,T?1ms;设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。

分析：f=200Hz时阻带衰减大于30dB，通过修改axis([0,fs/2,-80,10])为axis([200,fs/2,-1,1]) 发现通带波动rs满足<0.8。

16

bz =[0.0262 -0.1047 0.1570 -0.1047 0.0262] az =[1.0000 1.5289 1.6537 0.9452 0.2796] 系统函数为：

H(z)?0.0262-0.1047z1?1.5289z?1?1?0.1570z?2?2?0.1047z?3?3?0.0262z?4?4?1.6537z?0.9452z?0.2796z

（2）fc?0.2kHz，??1dB，fr?0.3kHz，At?25dB,T?1ms;分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。

17

bz1 =[0.0000 0.0002 0.0153 0.0995 0.1444 0.0611 0.0075 0.0002 0.0000 0] az1 =[1.0000 -1.9199 2.5324 -2.2053 1.3868 -0.6309 0.2045 -0.0450 0.0060 -0.0004] 因此脉冲响应不变法的系统函数为:

Himp(z)?0.0002z1?1.9199z?1?1?0.0153z?2?0.0995z?3?3?0.1444z?4?0.0611z?5?0.0075z?6?0.0002z?7?8

?0.0004z?9?2.5324z?2?2.2053z?1.3868z?4?0.6309z?5?0.2045z?6?0.0450z?7?0.0060zbz2 =[0.0179 0.1072 0.2681 0.3575 0.2681 0.1072 0.0179] az2 =[1.0000 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.0025] 因此双线性变换法的系统函数为:

Hbil(z)?0.0179?0.1072z1?0.6019z?1?1?0.2681z?2?2?0.3575z?3?3?0.2681z?4?4?0.1072z?5?5?0.0179z?6?6?0.9130z?0.2989z?0.1501z?0.0208z?0.0025z

分析：

脉冲响应不变法的N=9，双线性变换法的N=6，由图知它们都满足要求，但脉冲响应的衰减较快，双线性变换的过渡带窄一些，且阶数比脉冲小，容易实现。

实验四 FIR数字滤波器的设计

一、实验目的

1. 掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，

熟悉响应的计算机编程；

2. 熟悉线性相位FIR滤波器的幅频特性和相频特性； 3. 了解各种不同窗函数对滤波器性能的影响。 二、实验原理与方法

线性相位实系数FIR滤波器按其N值奇偶和h(n)的奇偶对称性分为四种： 1、h(n)为偶对称，N为奇数

18

H(ejω)的幅值关于ω=0，π，2π成偶对称。 2、h(n)为偶对称，N为偶数

H(ejω)的幅值关于ω=π成奇对称，不适合作高通。 3、h(n)为奇对称，N为奇数

H(ejω)的幅值关于ω=0，π，2π成奇对称，不适合作高通和低通。 4、h(n)为奇对称，N为偶数

H(ejω) ω=0、2π＝0，不适合作低通。 (一) 窗口法

窗函数法设计线性相位FIR滤波器步骤

确定数字滤波器的性能要求：临界频率{ωk}，滤波器单位脉冲响应长度N； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ejω)的幅频特性和相频特性；

求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ejω)按M(M远大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)； 选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应；

求H(ejω)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。

窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有： 矩形窗 w(n)=RN(n)； Hanning窗

；

Hamming窗 ；

Blackmen窗 ；

19

Kaiser窗 。 式中Io(x)为零阶贝塞尔函数。 （二）频率采样法

频率采样法是从频域出发，将给定的理想频率响应Hd(ejω)加以等间隔采样

然后以此Hd(k)作为实际FIR数字滤波器的频率特性的采样值H(k),即令

由H(k)通过IDFT可得有限长序列h(n)

将上式代入到Z变换中去可得

其中Φ(ω)是内插函数

（三）FIR滤波器的优化设计

FIR滤波器的优化设计是按照最大误差最小化准则，使所设计的频响与理想频响之间的最大误差，在通带和阻带范围均为最小，而且是等波动逼近的。 为了简化起见，在优化设计中一般将线性相位FIR滤波器的单位脉冲响应h(n)的对称中心置于n=0处，此时，线性相位因子α=0。当N为奇数，且N=2M+1，则

如希望逼近一个低通滤波器，这里M，ωp和ωs固定为某个值。在这种情况下有

定义一逼近误差函数：

20

本文来源：https://www.bwwdw.com/article/4yux.html