数字信号处理实验讲义2015-9-6

实验一 信号（模拟、数字）的输入输出实验

（常见离散信号产生和实现）

一、实验目的

1．加深对常用离散信号的理解；

2．掌握matlab中一些基本函数的建立方法。

二、实验原理

1.单位抽样序列

?(n)??

?1?0n?0n?0

在MATLAB中可以利用zeros()函数实现。

x?zeros(1,N);

x(1)?1;如果?(n)在时间轴上延迟了k个单位，得到?(n?k)即：

?1n?k ?(n?k)??

0n?0?2．单位阶跃序列

n?0?1 u(n)??

0n?0?在MATLAB中可以利用ones()函数实现。

x=ones(1,N)

3．正弦序列

x(n)?Asin(2?fn/Fs??)

在MATLAB中,

n=0:N-1;

x=A*sin(2*pi*f*n/Fs+phi)

1

4．复指数序列

x(n)?r?ej?n

在MATLAB中,

n=0:N-1;

x=r*exp(j*w*n) 5．指数序列

x(n)?an

在MATLAB中,

n=0:N-1;

x=a.^n

三、实验内容实现和图形生成

1．五种基本函数的生成 程序如下： (1)单位抽样序列 方法一：

% 单位抽样序列和延时的单位抽样序列 n=0:10;

x1=[1 zeros(1,10)];x2=[zeros(1,5) 1 zeros(1,5)]; subplot(1,2,1);

stem(n,x1);xlabel ('时间序列n');ylabel('振幅');title('单位抽样序列x1'); subplot(1,2,2);

stem(n,x2); xlabel('时间序列n');ylabel('振幅');title('延时了5的单位抽样序列');

2

方法二：

先在matlab中定义单位抽样序列： function [x,n]=dwxl(n1,n2,n0) n=[n1:n2]; x=(n==n0);

或：function [x,n]=impseq(n1,n2,n0) n=[n1:n2]; x=[(n-n0)==0]; 在运行命令： [x,n]=dwxl(-5,5,3)

stem(n,x);xlabel('n');title('单位抽样序列x'); 或：[x,n]=impseq(-5,5,3)

stem(n,x);xlabel('n');title('单位抽样序列x'); （2）单位阶跃序列 方法一： n=0:10; u=[ones(1,11)];

stem(n,u);xlabel ('时间序列n');ylabel('振幅');title('单位阶跃序列'); 所得的图形如下所示：

3

方法二;

先在matlab中定义单位阶跃序列： function [x,n]=jyxl(n1,n2,n0) n=[n1:n2]; x=(n>=n0);

或: function [x,n]=stepseq(n1,n2,n0)

n=[n1:n2]; x=[(n-n0)>=0];

在运行命令： [x,n]=jyxl(-5,5,3)

stem(n,x);xlabel('n');title('单位阶跃序列x'); 或：[x,n]=stepseq(-5,5,3)

stem(n,x);xlabel('n');title('单位阶跃序列x'); （3）正弦函数 n=1:30;

x=2*sin(pi*n/6+pi/3);

stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('正弦函数序列x=2*sin(pi*n/6+pi/3)');

(4)复指数序列

4

n=1:30; x=2*exp(j*3*n);

stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('复指数序列x=2*exp(j*3*n)'); 图形如下：

（5）指数序列 n=1:30; x=（-0.8）^n;

stem(n,x); xlabel ('时间序列n');ylabel('振幅');title('指数序列x=1.2.^n');

2、绘出信号x(n)?ezn，当s??(1/12?)j??、s?(1/12)?j时、z?112、66z?2?j?6、z?j?6时的信号实部和虚部图；

程序如下：

z1=-1/12+j*pi/6;z2=1/12+j*pi/6;z3=1/12;z4=2+j*pi/6;z5=j*pi/6; n=0:20;

x1=exp(z1*n);x2=exp(z2*n); x3=exp(z3*n);x4=exp(z4*n); x5=exp(z5*n); subplot(5,2,1);

stem(n,real(x1)); xlabel ('时间序列n');ylabel('实部');title('复指数z1=-1/12+j*pi/6时序列实部'); subplot(5,2,2);

stem(n,imag(x1)); xlabel ('时间序列n');ylabel('虚部');title('复指数z1=-1/12+j*pi/6时序列虚部');

5

subplot(5,2,3);

stem(n,real(x2)); xlabel ('时间序列n');ylabel('实部');title('复指数z2=1/12+j*pi/6时序列实部'); subplot(5,2,4);

stem(n,imag(x2)); xlabel ('时间序列n');ylabel('虚部');title('复指数z2=1/12+j*pi/6时序列虚部'); subplot(5,2,5);

stem(n,real(x3)); xlabel ('时间序列n');ylabel('实部');title('复指数z3=1/12时序列实部'); subplot(5,2,6);

stem(n,imag(x3)); xlabel ('时间序列n');ylabel('虚部');title('复指数z3=1/12时序列虚部'); subplot(5,2,7);

stem(n,real(x4)); xlabel ('时间序列n');ylabel('实部');title('复指数z4=2+j*pi/6时序列实部'); subplot(5,2,8);

stem(n,imag(x4)); xlabel ('时间序列n');ylabel('虚部');title('复指数z4=2+j*pi/6时序列虚部'); subplot(5,2,9);

stem(n,real(x5)); xlabel ('时间序列n');ylabel('实部');title('复指数z5=j*pi/6时序列实部'); subplot(5,2,10);

stem(n,imag(x5)); xlabel ('时间序列n');ylabel('虚部');title('复指数z5=j*pi/6时序列虚部');

6

由上图的实部部分可以看出，Z=pi/6时，序列周期为12。计算序列周期为2*6＝12。实验和理论相符。

3．绘出信号x(n)?1.5sin(2?*0.1n)的频率是多少？周期是多少？产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期? 程序如下： n=0:40;

x1=1.5*sin(2*pi*0.1*n);x2=sin(0.9*n); subplot(1,2,1);

stem(n,x1); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x1=1.5*sin(2*pi*0.1*n)'); subplot(1,2,2);

stem(n,x2); xlabel ('时间序列n');ylabel('振幅');title('正弦序列x2=sin(0.9*n)'); 运行结果如下：

7

由上图看出：x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。理论计算中对第一个，N＝2*ｐｉ/（0.1*ｐｉ）＝10，第二个0.9不是pi的倍数，所以不是周期的。因此可以看出，实验结果和理论相符。

4．使用帮助功能学习 square(方波), sawtooth(锯齿波)和sinc 函数，并绘图。 （1）方波绘图程序如下： %用square t=-2*pi:0.001:2*pi; x=square(t); plot(t,x);

xlabel('t'),ylabel(' x=square(t)');

(2)三角波绘图程序如下： %用Sawtooth t=-2*pi:0.001:2*pi; y=sawtooth(t); plot(t,y);

xlabel('t'),ylabel(' y=sawtooth(t);');

8

（3）sinc函数绘图程序如下： t=-pi:0.001:pi; x=sinc(t); plot(t,x);

xlabel('t'),ylabel('sinc(t);');

四、问题讨论与总结：

1．离散正弦序列的性质：

离散正弦序列就是一个连续的正弦信号被一系列冲激函数采样后的结果，原连续正弦函数一定是周期的，但采样后的离散序列却不一定是周期的。对于离散序列x=sin(n*w)来说， 只有当2*pi/ｗ是一个有理数时，也就是说当ｗ是pi的倍数时，此离散序列才是周期的。所以在本实验中x1=1.5*sin(2*pi*0.1*n)的周期是10，而x2=sin(0.9*n)是非周期的。因为0.9不是pi的倍数。

2．离散复指数序列性质：

对于离散复指数函数ｘ=a*exp(z*ｎ)，只有当z是纯虚数，且纯虚数的系数是pi的倍数时，才是周期的。其它情况下均不是。这个性质由本次实验中的五个函数的图像可以被证明。

9

实验二 离散系统时域分析

一、实验目的

1．学习MATLAB语言的编程和调试技巧； 2．差分方程的求解；

3．掌握笔算离散卷积方法和MATLAB语言实现。

二、 实验内容

时域中，离散时间系统对输入信号或延迟信号进行运算处理，生成具有所需特性的输出信号。本实验通过MATLAB仿真一些简单的离散时间信号和系统，并研究其时域特性。涉及到离散时间信号、离散时间系统、系统性质及线性卷积等知识点。

三、实验原理与方法和手段

一个离散时间系统，输入信号为x(n)，输出信号为y(n)，运算关系用T[﹒]表示，则输

入与输出的关系可表示为y(n)＝T[x(n)]。

在《信号与系统》和《数字信号处理》课程中，我们知道描述线性移不变离散时间系统的数学模型是常系数差分方程，它与系统的结构流图之间可以互相推导。迭代解法（也称递推解法）是求解差分方程的最简单也最适用的方法，也是实现数字滤波器的一种基本方法。 差分方程通式为：

?ay?n?k???bx?n?m?

kmk?0m?0NM?a[k]y(n?k)??b[r]x(n?r)k?0r?0NNx(n)与y(n)分别为系统的激励和响应。

MATLAB以函数filter(b, a , x)，来计算在给定输入和差分方程系数时求差分方程的数值解。b,a分别为系统方程的系数向量。x是输入序列。

线性时不变系统的输入输出关系可通过单位脉冲响应h(n)表示： y(n)?x?n?*h(n)?式中*表示卷积运算。

可物理实现的线性时不变系统是稳定的、因果的。这种系统的单位脉冲响应是因果的

m????x(m)h(n?m)

?? 10

（单边）且绝对可和的，即：h(n)=0,n<0；

n????h(n)??。在MATLAB语言中采用conv

??实现卷积运算即：y=conv(x,h)，它默认从n=0开始。

四、实验组织运行要求

1．学生在进行实验前必须进行充分的预习，熟悉实验内容； 2．学生根据实验要求，读懂并理解相应的程序；

3．学生严格遵守实验室的各项规章制度，注意人身和设备安全，配合和服从实验室人员管理；

4．教师在学生实验过程中予以必要的辅导，独立完成实验； 5．采用集中授课形式。

五、实验条件

1．具有WINDOWS 98/2000/NT/XP操作系统的计算机一台； 2.。MATLAB编程软件。

六、实验步骤

在“开始--程序”菜单中，找到MATLAB程序，运行启动；进入MATLAB后 ，首先熟悉界面；在Command Window中输入参考程序，并执行；记录运行结果图形，并与笔算结果对照。（MATLAB的使用请参考附录） 具体步骤如下：

1．已知某一系统方程为： y(n)+3y(n-1)+2y(n-2)=x(n)

（1）利用MATLAB计算系统的脉冲响应h(n)并画出波形图,其中n=(-2,10); （2）若x(n)=u(n)，利用MATLAB计算系统的零状态响应y(n)并画出波形图。 2．设某LSI的单位脉冲响应和h(n)=()u(n), （1）判断此系统是否可实现；

（2）当输入为x(n)?3[u(n)?u(n?5)]时，利用卷积和手工计算此LSI的输出y(n)； （3）用MATLAB实现第二步，并画出图形(h(n)的n取值范围自定)。

3．x(n)=[2,3,1,1,2,-1,0,3]，-2≤n≤5；h(n)=[2,4,1,-2,0,-1]，-3≤n≤2，手工计算和MATLAB

11

n12n

计算卷积y(n)=x(n)*h(n)。

七、思考题

结合《信号与系统》课程所学，思考离散时间系统的线性卷积公式与连续时间系统的卷积公式的异同？

八、实验报告要求

1．报告中要给出实验的MATLAB程序； 2．简述实验目的和原理；

3．给出用笔算时差分方程解、卷积和conv计算线性卷积对照图； 4．给出收获和体会。

九、参考程序

1．程序1

系统微分方程y(n)-y(n-1)+0.9y(n-2)=x(n)，计算并画出系统单位冲激响应h(n),n=(-5,10)。 %xh08 n=[-5:10]; b=[1]; a=[1 -1 0.9]; x=impseq(-5,10,0); h=filter(b,a,x); stem(n,h, '.');

xlabel('时间序号n'); ylabel('h(n)'); title('单位冲激响应');grid on 2．程序2

设某LSI的单位脉冲响应h(n)?0.8u(n) （1）判断此系统是否可实现；

（2）当输入为矩形脉冲x(n)?2[u(n)?u(n?10)]时，手工计算此LSI的输出y(n)； %xh09 n1=0:9;

12

nn

x=2.^n1; N1=length(x); n2=0:19;

h=0.8.^n2; N2=length(h); y=conv(x,h); N=N1+N2-1; n=0:N-1;

subplot(2,2,1);stem(n1,x,'.');xlabel('时间序号n');title('输入序列x(n)'); grid on; subplot(2,2,2); stem(n2,h,'.');xlabel('时间序号n');title('单位取样序列h(n)');grid on; subplot(2,2,3);stem(n,y,'.');xlabel('时间序号n');title('输出序列y(n)');grid on; 3．程序3

如果x(n)、h(n)的起点不为0，则采用conv_m计算卷积； 编写conv_m函数：

function[y,ny]=conv_m(x,nx,h,nh) %改进卷积程序 nyb=nx(1)+nh(1);

nye=nx(length(x))+nh(length(h)); ny=[nyb:nye]; y=conv(x,h); 在命令窗口输入： x=[3,11,7,0,-1,4,2];nx=[-3:3]; h=[2,3,0,-5,2,1];nh=[-1:4]; [y,ny]=conv_m(x,nx,h,nh)

stem(ny,y,'.');xlabel('时间序号n');title('卷积和y(n)=x(n)*h(n)');grid on;

4．function [x,n]=impseq(n1,n2,n0) function [x,n]=stepseq(n1,n2,n0) n=[n1:n2]; n=[n1:n2];

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

13

实验三 FFT频谱分析及应用

一、实验目的

1．通过实验加深对FFT的理解；

2．熟悉应用FFT对典型信号进行频谱分析的方法。

二、实验内容

使用MATLAB程序实现信号频域特性的分析。涉及到离散傅立叶变换(DFT)、快速傅立叶变换(FFT)及信号频率分辨率等知识点。

三、实验原理与方法和手段

在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换

(DFT)进行分析。DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。

有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。

在MATLAB信号处理工具箱中的函数fft（x,N）,可以用来实现序列的N点快速傅立叶变换。经函数fft求得的序列一般是复序列，通常要求出其幅值和相位。MATLAB中提供了求复数的幅值和相位的函数：abs、angle，这些函数一般和fft同时使用。

四、实验组织运行要求

1．学生在进行实验前必须进行充分的预习，熟悉实验内容； 2．学生根据实验要求，读懂并理解相应的程序；

3．学生严格遵守实验室的各项规章制度，注意人身和设备安全，配合和服从实验室人员管理；

4．教师在学生实验过程中予以必要的辅导，独立完成实验； 5．采用集中授课形式。

五、实验条件

1．具有WINDOWS 98/2000/NT/XP操作系统的计算机一台；

14

2.。MATLAB编程软件。

六、实验步骤

在“开始--程序”菜单中，找到MATLAB程序，运行启动；

进入MATLAB后 ，在Command Window中输入实验程序，并执行； 记录运行结果图形，作分析。 具体步骤如下：

1． 用FFT 进行典型信号的频谱分析 高斯序列：

??

改变参数p、q,分析参数的变化对频谱的影响。

2．模拟信号x(t)?2sin(4?t)?5cos(8?t)，以时间以0.01s为间隔即t?0.01n进行采样，其中n=(0:N-1)，求：

（1）N＝40点FFT的幅度频谱，从图中能否观察出信号的2个频谱分量?

（2）提高采样点数，如N＝128及1024，再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的2个模拟频率和数字频率各为多少？FFT频谱分析结果与理论上是否一致？

3．有限长序列x(n)={2,1,0,1,3};h(n)={1,3,2,1}，试利用FFT实现由DFT计算线性卷积，并与线性卷积直接计算（conv）的结果进行比较。

七、实验报告要求

1．报告中要给出实验的MATLAB程序； 2．简述实验目的和原理；

3．按实验步骤附上实验信号序列和幅频特性曲线，分析所得到的图形，回答每一步提出的问题。

八、参考程序

1．程序1：

15

%xh10

?T进行高斯序列的频谱分析

n=0:15; %定义序列的长度是15 p=8; q=2; x=exp(-1*(n-p).^2/q); %利用fft 函数实现付氏变换 close all; subplot(3,2,1); stem(x,’.’);subplot(3,2,2);stem(abs(fft(x)),’.’);grid on; p=8; q=4; x=exp(-1*(n-p).^2/q); %改变信号参数，重新计算 subplot(3,2,3);stem(x,’.’);subplot(3,2,4); stem(abs(fft(x)),’.’);grid on; p=8; q=8; x=exp(-1*(n-p).^2/q);

subplot(3,2,5);stem(x,’.’);subplot(3,2,6); stem(abs(fft(x)),’.’);grid on; 2．程序2： %xh11 N=40;n=0:N-1; t=0.01*n;

x=2*sin(4*pi*t)+5*cos(8*pi*t); k=0:N/2;w=2*pi/N*k; X=fft(x,N);

magX=abs(X(1:N/2+1));

subplot(2,1,1);stem(n,x,'.');title('signal x(n)'); subplot(2,1,2);plot(w/pi,magX);title('FFT N=40'); xlabel('w (unit :pi)');ylabel('|X|');grid 3．程序3： %xh12

% 用FFT实现由DFT计算线性卷积 x=[1 2 0 1];h=[2 2 1 1]; L=length(x)+length(h)-1; XE=fft(x,L);HE=fft(h,L); y1=ifft(XE.*HE);

%画出由圆周卷积计算线性卷积结果及误差 k=0:L-1;

subplot(221);stem(k,real(y1),’.’);axis([0 6 0 7]);

16

title('圆周卷积');

xlabel('Time index k');ylabel('Amplitude');grid on; y2=conv(x,h);

subplot(222);stem(k, y2,’.’);axis([0 6 0 7]); title('线性卷积');

xlabel('Time index k');ylabel('Amplitude');grid on; error=y1-y2; subplot(223); stem(k,abs(error));

xlabel('Time index k');ylabel('Amplitude'); title('Error Magnitude');grid on;

17

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

一、实验目的

1．掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理和方法；

2．观察双线性变换法和脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别。

二、实验内容

使用MATLAB编写程序，实现IIR数字滤波器的设计。涉及脉冲响应不变法和双线性变换法设计IIR数字滤波器的方法、不同设计方法得到的IIR滤波器频域特性异同等知识点。

三、实验原理与方法和手段

1．脉冲响应不变法

所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应ha(t)的采样值，即：h(n)?ha(t)t?nT?ha(nt)，其中，T为采样周期。

在脉冲响应不变法中，模拟角频率和数字角频率的变换关系为：???T，可见，?和ω之间的变换关系为线性的。

在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：

[B,A]=impinvar(b,a,fs1) [B,A]=impinvar(b,a)

其中，b、a分别为模拟滤波器的分子和分母多项式系数向量；fs1为采样频率（Hz），缺省值fs=1Hz；B、A分别为数字滤波器分子和分母多项式系数向量。 2．双线性变换法：

21?z?1由于s平面和z平面的单值双线性映射关系为s＝，其中T为采样周期。因此，?1T1?z若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。

在双线性变换中，模拟角频率和数字角频率的变换关系为：??ω之间的变换关系为非线性的。

在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，

18

2?tg，可见，Ω和T2

调用格式为：[B,A]=bilinear(b,a,fs1) 3．数字滤波器设计

（1）定技术指标转换为模拟滤波器设计性能指标。

（2）估计满足性能指标的模拟相应滤波器性能阶数和截止频率。 利用MATLAB中buttord、cheb1ord、cheb2ord、ellipord等函数， 调用格式如：[N,?c]?buttord(?p,?s,?p,?s,'s')

其中，Ωp为通带边界频率，rad/s；Ωs为阻带边界频率，rad/s;αp为带通波动，dB;αs为阻带衰减，dB；‘s’表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Ωc为模拟滤波器的截止频率(-3dB频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。 （3）设计模拟滤波器。

MATLAB信号处理工具箱提供了模拟滤波器设计的完全工具函数：butter、cheby1，cheby2、ellip、besself。用户只需一次调用就可完成低通、高通、带通、带阻滤波器设计。 调用格式如：[b,a]=butter(N,Ωc,’ftype’,‘s’)，其中，’ftype’为滤波器类型： ‘high’表示高通滤波器，截止频率为Ωc；

‘stop’表示带阻滤波器，Ωc=[Ω1 Ω2] (Ω1<Ω2)； ‘ftype’缺省时表示为低通或带通滤波器。

（４）利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。

四、实验组织运行要求

1．学生在进行实验前必须进行充分的预习，熟悉实验内容； 2．学生根据实验要求，编写相应的程序；

3．学生严格遵守实验室的各项规章制度，注意人身和设备安全，配合和服从实验室人员管理；

4．教师在学生实验过程中予以必要的辅导，独立完成实验； 5．采用集中授课形式。

五、实验条件

1．具有WINDOWS 98/2000/NT/XP操作系统的计算机一台； 2．MATLAB编程软件。

19

六、实验步骤

在“开始--程序”菜单中，找到MATLAB程序，运行启动；

进入MATLAB后 ，在Command Window中输入自己编写的主程序，并执行；记录运行结果图形，作分析对比。 具体步骤如下：

1、查看帮助文件，了解相关函数的调用格式。

2、设计模拟巴特沃斯低通滤波器，fp=5kHz，αp=2dB，fs=12kHz，αs=30dB。

3、设计模拟巴特沃斯高通滤波器，fp=200Hz，fs=100Hz，幅度特性单调下降，fp处最大衰减为3dB，阻带最小衰减αs=15dB。

4、fp=0.2kHZ，αp=1dB，fs=0.4kHZ，αs=30dB，采样间隔T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器

(1)观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量； (2)比较两种方法的优缺点；

5、设计一数字高通滤波器，它的通带为700～1000Hz，通带内容许有2dB的波动，阻带内衰减在小于500Hz的频带内至少为30dB，采样频率为2000Hz。

6、利用MATLAB编程设计一个Butterworth数字带通滤波器，指标要求如下： 通带边缘频率：ωp1=0.5?rad，ωp2=0.7?rad，通带峰值起伏：αp≤2dB; 阻带边缘频率：ωs1=0.35?rad，ωs2=0.95?rad，最小阻带衰减：αs≥35dB。 分别用脉冲响应不变法和双线性变换法进行IIR数字滤波器的设计。T=0.2ms

七、实验报告要求

1．报告中要给出实验的MATLAB程序，了解每个语句作用； 2．简述实验目的和原理；

3．按实验步骤附上所设计滤波器的H(z)及相应的幅频特性曲线定性分析得到的图形，判断设计是否满足要求；

4．总结双线性变换法和脉冲不变法的特点和区别； 5．收获与建议。

八、MATLAB示例

例1：设计模拟巴特沃斯低通滤波器，fp=300Hz，αp=1dB，fs=800Hz，αs=20dB。

20

%xh13模拟低通滤波器技术指标 αp=1; αs=20; fp=300; fs=800; Ωp=2*pi*fp; Ωs=2*pi*fs;

%设计模拟巴特沃斯低通滤波器

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s')

[b,a]=butter(N,Ωc,'s') %设计模拟巴特沃斯低通滤波器，Ωp为通带边界频率，rad/s；Ωs为阻带边界频率，rad/s;αp为通带最大衰减，dB；αs为阻带最小衰减，dB；'s'表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Ωc为模拟滤波器的截止频率(-3dB频率)，rad/s；b、a分别为模拟滤波器的系统函数分子和分母多项式系数向量； [H,Ω]=freqs(b,a); %求模拟滤波器的频率响应 %绘制频响幅度谱

plot(Ω/2/pi,20*log10(abs(H))); %横轴为频率，单位：HZ;纵轴频响幅度，单位：dB axis([0,1500,-50,0]);

xlabel('频率Hz');ylabel(' H幅值dB');

例2：设计模拟巴特沃斯高通滤波器，fp=800Hz，αp=1dB，fs=300Hz，αs=20dB。 %xh14模拟高通滤波器技术指标 αp=1; αs=20; fp=800; fs=300; Ωp=2*pi*fp; Ωs=2*pi*fs;

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s')

[b,a]=butter(N,Ωc,'high','s') %设计模拟巴特沃斯高通滤波器，Ωp为通带边界频率，rad/s；Ωs为阻带边界频率，rad/s；αp为通带最大衰减，dB；αs为阻带最小衰减，dB；'s'表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Ωc为模拟滤波器的

21

截止频率(-3dB频率)，rad/s；b、a分别为模拟高通滤波器的系统函数分子和分母多项式系数向量；

[H,Ω]=freqs(b,a);

plot(Ω/2/pi,20*log10(abs(H))); axis([0,1000,-50,0]);

xlabel('频率Hz');ylabel(' H幅值dB');

例3：fp=0.１kHZ,αp=1dB，fs=0.3kHZ，αs=25dB，T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器。 %xh15采用冲击响应不变法 αp=1; αs=25; fp=100; fs=300; Ωp=2*pi*fp;

Ωs=2*pi*fs; %数字滤波器技术指标要求转化成模拟滤波器技术指标要求 fs1=1000; %采样频率 %设计模拟滤波器

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s'); [b,a]=butter(N,Ωc,'s');

[B,A]=impinvar(b,a,fs1) %用冲击响应不变法将模拟滤波器变换成数字滤波器。B、A分别为数字滤波器的系统函数分子和分母多项式系数向量； [H1,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应 %绘制数字滤波器频响幅度谱

subplot(211);plot(w*fs1/2/pi,20*log10(abs(H1))); axis([0,1000,-100,0]);

xlabel('频率Hz');ylabel(' H1幅值dB'); %采用双线性变换法 %数字滤波器的技术指标要求

wp=2*pi*fp/fs1; %数字滤波器通带边界频率 ws=2*pi*fs/fs1; %数字滤波器阻带边界频率

22

%变换为同类型模拟滤波器的技术指标要求

Ωp=2*fs1*tan(wp/2); %同类模拟滤波器通带边界频率 Ωs=2*fs1*tan(ws/2); %同类模拟滤波器阻带边界频率 %设计同类型模拟滤波器

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s'); [b,a]=butter(N,Ωc,'s');

%用双线性变换法将模拟滤波器变换成数字滤波器 [B,A]=bilinear(b,a,fs1)

[H2,w]=freqz(B,A,'whole'); %求数字滤波器的频率响应 %绘制数字滤波器频响幅度谱

subplot(212);plot(w*fs1/2/pi,20*log10(abs(H2))); axis([0,1000,-100,0]);

xlabel('频率Hz');ylabel('H2幅值dB');

例4：设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。 %xh16

Ωp=2*1000*tan(2*pi*400/(2*1000)); Ωs=2*1000*tan(2*pi*317/(2*1000)); [N,Ωc]=cheb1ord(Ωp,Ωs,0.5,19,'s'); [b,a]=cheby1(N,0.5,Ωc,'high','s'); [B,A]=bilinear(b,a,1000); [H,w]=freqz(B,A);

f=w/pi*500; %f等于数字角频率除以2?Ts plot(f,20*log10(abs(H))); axis([0,500,-80,10]); grid;

xlabel('频率Hz'); ylabel('幅度/dB');

例5： 利用MATLAB编程设计一个Butterworth数字带通滤波器，指标要求如下： 通带边缘频率：ωp1=0.45?rad，ωp2=0.65?rad，通带峰值起伏:αp≤1dB; 阻带边缘频率：

23

ωs1=0.3?rad，ωs2=0.8?rad，最小阻带衰减：αs≥40dB。分别用脉冲响应不变法和双线性变换法进行IIR数字滤波器的设计。Ts=0.125ms 1.双线性变换法（巴特沃兹原型）： %xh17 fs1=8000;

Ωs1=2*fs1*tan(0.3*pi/2); Ωs2=2*fs1*tan(0.8*pi/2); Ωp1=2*fs1*tan(0.45*pi/2); Ωp2=2*fs1*tan(0.65*pi/2); Ωs=[Ωs1 Ωs2]; Ωp=[Ωp1 Ωp2]; αp=1; αs=40;

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s'); [b,a]=butter(N,Ωc,'s'); [B,A]=bilinear(b,a,fs1); [H,w]=freqz(B,A); f=w/2/pi*fs1; subplot (2,1,1);

plot(f,20*log10(abs(H))); axis([0,4000,-60,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB'); subplot(2,1,2); plot(f,angle(H));

grid; xlabel('频率/Hz') ;ylabel('相位'); 2.脉冲响应不变法（巴特沃兹原型）： %xh18 Fs1=8000;

Ωs1=0.3*pi*fs1; Ωs2=0.8*pi*fs1; Ωp1=0.45*pi*fs1; Ωp2=0.65*pi*fs1; Ωs=[Ωs1 Ωs2]; Ωp=[Ωp1 Ωp2]; αp=1; αs=40;

24

[N,Ωc]=buttord(Ωp,Ωs,αp,αs,'s'); [b,a]=butter(N,Ωc,'s'); [B,A]=impinvar(b,a,fs1); [H,w]=freqz(B,A); f=w/pi*4000; subplot(2,1,1);

plot(f,20*log10(abs(H))); axis([0,4000,-80,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB'); subplot(2,1,2); plot(f,angle(H));

grid; xlabel('频率/Hz') ;ylabel('相位');

25

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

一、实验目的

1．掌握用窗函数法设计FIR数字滤波器的原理和方法； 2．熟悉线性相位FIR滤波器的幅频特性和相频特性； 3．了解不同窗函数对滤波器性能的影响。

二、实验内容

使用MATLAB编写程序，实现FIR数字滤波器的设计。涉及窗函数法和频率采样法设计FIR数字滤波器的方法、线性相位FIR滤波器的幅频特性和相频特性的特点、窗函数选择及其对滤波器性能的影响等知识点。

三、实验原理与方法和手段

窗函数设计法是一种把一个长序列变成有限长的短序列的设计方法，是在时域进行的。用窗函数法设计FIR数字滤波器时，先根据W c和N 求出相应的的理想滤波器的单位脉冲响应hd(n)。

1hd(n)?2??wc?wcHd(ejw)ejwndw?sin(wc(n?a))?(n?a)

因为hd(n)一般是非因果的，且无限长，物理上是不可实现的。为此可选择适当的窗函数w(n)截取有限长的hd(n)， 即h(n)= hd(n)w(n)，只要阶数足够长，截取的方法合理，总能够满足频域的要求。

实际中常用的窗函数有矩形（Boxcar）窗、三角（Bartlett）窗、汉宁(Hanning)窗、汉明(Hamming) 窗和布莱克曼（Blackman）窗。这些窗函数各有优缺点,所以要根据实际情况合理选择窗函数类型。

1．窗函数法设计线性相位FIR滤波器的一般步骤为： （1）确定理想滤波器Hd(e（2）由Hd(e

j?j?)的特性；

)求出hd(n)；

26

（3）根据过渡带宽度和阻带最小衰减，借助窗函数确定窗的形式及N的大小，即选择适当的窗函数，并根据线性相位条件确定窗函数的长度N；在MATLAB中，可由w=boxcar(N)（矩形窗）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗）等函数来实现窗函数设计法中所需的窗函数。 （4）由h(n)?hd(n)?w(n), 0≤n≤ N-1，得出单位脉冲响应h(n); （5）对h(n)作离散时间傅立叶变换，得到H(ej?)。

2．在MATLAB中，可以用b=fir1(N,Wn,‘ftype’,taper) 等函数辅助设计FIR数字滤波器。N代表滤波器阶数；Wn代表滤波器的截止频率(归化频率)，当设计带通和带阻滤波器时，Wn为双元素相量；ftype代表滤波器类型，如‘high’高通，‘stop’带阻等；taper为窗函数，默认为海明窗，窗函数实现需要用窗函数blackman, hamming,hanning chebwin, kaiser产生。

四、实验组织运行要求

1．学生在进行实验前必须进行充分的预习，熟悉实验内容； 2．学生根据实验要求，编写相应的程序；

3．学生严格遵守实验室的各项规章制度，注意人身和设备安全，配合和服从实验室人员管理；

4．教师在学生实验过程中予以必要的辅导，独立完成实验； 5．采用集中授课形式。

五、实验条件

具有WINDOWS 2000/XP操作系统的计算机一台，安装MATLAB软件

六、实验步骤

在“开始--程序”菜单中，找到MATLAB程序，运行启动；

进入MATLAB后 ，在Command Window中输入自己编写的主程序，并执行； 记录运行结果图形，作分析对比。

具体步骤如下：

27

1．用窗函数法设计一线性相位FIR低通滤波器，要求通带截止频率wc=

?， 4（1）选择一个合适的窗函数（如hamming窗），取单位冲击响应h（n）的长度N＝15，观察所设计滤波器的幅频特性，记录阻带最小衰减；

（2）取N＝45，重复上述设计，观察幅频和相频特性的变化，分析长度N变化的影响； （3）保持N＝45不变，改变窗函数（如hamming窗变为blackman窗），观察并记录窗函数对滤波器幅频特性的影响。

2．针对一个含有35Hz、50Hz和70Hz的混和正弦波信号，设计一个FIR带通滤波器。 参数要求:采样频率fs=400Hz，通带下限截止频率fc1=40Hz，通带上限截止频率fc2=60Hz，过渡带宽8Hz，通阻带波动0.01，采用凯塞窗设计。

七、实验报告要求

1.报告中要给出实验的MATLAB程序，了解语句作用； 2.简述实验目的和原理；

3.按实验步骤附上所设计滤波器的h(n)及相应的幅频和相频特性曲线，比较它们的性能，说明不同的窗函数对滤波器性能的影响； 4.收获和建议。

八、MATLAB示例

1.用窗函数法设计一线性相位FIR低通滤波器，要求通带截止频率wc=（n）的长度N＝21。绘出h(n)及其幅频响应曲线.

?，单位冲击响应h3sin((n?a))3 设计分析:理想滤波器单位冲击响应为hd(n)?；为了满足线性相位FIR滤?(n?a)N?1波器条件h(n)= h(N-1-n)，要求a==10。

2%xh19

N=21;wc=pi/3; %理想低通滤波器参数 n=0:N-1; a=(N-1)/2;

hdn=sin(wc*(n-a))/pi./(n-a); %计算理想低通滤波器单位冲击响应hd(n) if rem(N,2)~=0 hdn(a+1)=wc/pi; end %N为奇数时,处理n=a点的0/0型

28

?

wn=hamming(N); %hamming窗 hn=hdn.*wn';

subplot(121); %绘出h(n)及幅频特性曲线 stem(n,hn,'.');

xlabel('n');ylabel('h(n)'); title('hamming窗设计的h(n)'); grid on; subplot(122); hw=fft(hn,512); w=2*[0:511]/512;

plot(w,20*log10(abs(hw)));

xlabel('w/pi');ylabel('Magnitude(dB)'); title('hamming窗设计的幅频特性'); grid on;

2．针对一个含有5Hz、15Hz和30Hz的混和正弦波信号，设计一个FIR带通滤波器。 参数要求:采样频率fs=100Hz，通带下限截止频率fc1=10Hz，通带上限截止频率fc2=20Hz，过渡带宽6Hz，通阻带波动0.01，采用凯塞窗设计。 %xh20

fs=100; T=1/fs;

[n,wc,beta,ftype]=kaiserord([7 13 17 23],[0 1 0], [0.01 0.01 0.01],fs) wn=kaiser(n+1,beta); %使用kaiser窗函数，(n+1)为窗宽 b=fir1(n,wc,wn); %使用标准频率响应的加窗设计函数fir1 freqz(b,1,512); %数字滤波器频率响应，a=1 t=0:T:1;

s=sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30); sf=filter(b,1,s); %对信号s进行滤波 figure

subplot(2,1,1); plot(t,s);grid on; subplot(2,1,2); plot(t,sf);grid on

针对一个含有15Hz和40Hz的混和正弦波信号，设计一个FIR高通滤波器。

29

%高通滤波器 fs=200;

[n,Wn,beta,ftype]=kaiserord([27 33],[0 1], [0.01 0.01],200) window=kaiser(n+1,beta); %使用kaiser窗函数，(n+1)为窗宽 b=fir1(n,Wn,'high',window); %使用标准频率响应的加窗设计函数fir1 freqz(b,1,512); %数字滤波器频率响应 t=(0:100)/fs;

s=sin(2*pi*t*15)+sin(2*pi*t*40);

sf=filter(b,1,s); %对信号s进行滤波 figure

subplot(2,1,1); plot(t,s);grid on; subplot(2,1,2); plot(t,sf);grid on

30

实验六 数字信号处理综合实验

一、实验目的

1．学习MATLAB的离散仿真环境和用MATLAB建立离散系统的方法；; 2．通过实验操作熟悉滤波器结构； 3．设计一数字滤波器。

二、实验内容

1．用MATLAB仿真离散系统；用直接Ⅰ型、直接Ⅱ型结构实现滤波器； 2．设计一个综合数字滤波器。

三、实验原理与方法和手段

FIR滤波器、IIR滤波器的各种结构；MATLAB中提供的模块组建各种滤波器。

用MATLAB仿真离散系统：已知某系统的差分方程为y(n)=x(n)+2.5y(n-1)-y(n-2), 其冲激响应h(n)=[4/3(2)n-1/3(0.5)n]u(n),用MATLAB仿真差分方程，求其冲激响应，并把结果与h(n)比较。 具体操作步骤：

1．在MATLAB命令窗口中输入simulink并回车以打开仿真模块库；

2．用单位冲激信号、加法器、延时器、数乘器、输出接口等模块按下图连接； 3．点击仿真窗口工具条中的图标

开始仿真；

4．回到MATLAB命令窗口中，输入以下程序，比较结果,同时在MATLAB workspace窗口得到y的数值解。 %xh14 n=0:10;

h=(4/3)*2.^n-(1/3)*(0.5).^n; subplot(2,2,1);

31

stem(n,h);%作出理论上的冲激响应 subplot(2,2,2);

stem(n,yout1);%作出实验中直接型的冲激响应 subplot(2,2,3);

stem(n,yout2);%作出实验中级联型的冲激响应 subplot(2,2,4);

stem(n,yout3);%作出实验中并联型的冲激响应

四、实验组织运行要求

1．学生在进行实验前必须进行充分的预习，熟悉实验内容； 2．学生根据实验要求，编写相应的程序；

32

3．学生严格遵守实验室的各项规章制度，注意人身和设备安全，配合和服从实验室人员管理；

4．教师在学生实验过程中予以必要的辅导，独立完成实验； 5．采用集中授课形式。

五、实验条件

具有WINDOWS 2000/XP操作系统的计算机一台，安装MATLAB软件

六、实验步骤

1．用直接Ⅰ型结构实现滤波器： 已知某滤波器的差分方程：

y(t)?x(n)?响应。

131x(n?1)?y(n?1)?y(n?2),用直接Ⅰ型结构实现该滤波器，求冲激348（1）在MATLAB file中新建一个model文件；

（2）用单位冲激信号、加法器、延时器、数乘器、输出接口等模块,按下图连接各模块；

（3）点击仿真窗口工具条中的图标

开始仿真；

（4）回到MATLAB命令窗口中，输入以下程序，比较结果,同时在MATLAB workspace窗口得到y的数值解。

33

n=0:20;

stem(n,yout1);

2．用直接Ⅱ型结构实现滤波器：

用直接Ⅱ型结构实现步骤1中的滤波器，求冲激响应。按下图新建一个model文件，具体步骤同上。

n=0:20;

subplot(211);stem(n,yout1);%直接Ⅰ型冲激响应 subplot(212);stem(n,yout2); %直接Ⅱ型冲激响应

七、实验报告要求

1．简述实验目的和原理；

2．用级联型和并联型实现滤波器，必须给出仿真模块图和冲激响应图及数值； 3．收获和建议。

34

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