数字信号处理作业

1 利用MATLAB设计数字滤波器

随着计算机和信息科学的极大发展，信号处理己经逐步发展为

一门独立的学科，是信息科学的重要组成部分，在语音处理、图像处理、雷达、航空航天、地质勘探、通信、生物医学等众多领域得到了广泛的应用。信号是数字信号处理领域最基本最重要的概念，信号是信息的载体是信息的物理体现。而数字滤波器作为信号处理一项关键技术是数字信号处理的重要基础，在对信号的过滤、检测、与参数估算等处理过程中，它是使用最为广泛的一种线性系统。数字滤波器按照其冲激响应函数的时域特性，可分为无限长冲激响应(infinite impulse response, IIR)滤波器和有限长冲激响应(finite impulse Response, FIR)滤波器。在满足相同指标下，IIR滤波器的阶数明显小于FIR，硬件实现容易目大大减少了运算量，在不要求严格线性相位的情况下，IIR滤波器的应用相当广泛。

Matlab是一种交互式的以矩阵为基础的软件，它用于科学与工程项日的计算与可视化，它只需要其它编程语言的几分之一的时间即可以解决复杂的数值计算问题。它的强有力也在于那些相对简单的编程功能和提供的非常方便简单的不同学科的工具箱。在设计数字滤波器的也相当方便。

1

2 IIR数字滤波器的设计

相对于FIR数字滤波器，IIR数字滤波器最突出的有点是：IIR数字滤波器能够以更低的阶数n满足相同的技术参数要求。

2.1 IIR数字滤波器的基本概念

若N阶递归型数字滤波器的差分方程为

y(n)bix(nr)aiy(ni)

i0i1MN则IIR滤波器的系统函数

H(Z)bzrr1Nk1Mr

1akzk 从以上的系统函数可知，设IIR滤波器的任务就是通过计算寻求一个因果、物理上可实现的系统函数H(Z)，使其频率响应

H(jw)满足所希望得到的频域指标，即符合给定的通带截止频率、

阻带截止通带衰减和阻带衰减。不难看出，数字滤波器与模拟滤波器的设计思路相仿，其设计实质也是寻找一组系数{a,b},去逼近所要求的频率响应，使其在性能上满足预定的技术要求；不同的是模拟滤波器的设计是在S平面上用数学逼近法去寻找近似的所需特性H(S)，而数字滤波器则是在Z平面寻找合适的H(z)。IIR数字滤波器的单位响应是无限长的，而模拟滤波器一般都具有无限长的单位脉冲响应，因此与模拟滤波器相匹配。由于模拟滤波器的设计在理论上已十分成熟，因此数字滤波器设计的关键是将

2

H(S)→H(Z)，即利用复值映射将模拟滤波器离散化。

2.2 双线性变换法设计数字滤波器

基本设计过程：

（1）

将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；

（2） （3）

设计过度模拟滤波器；

将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。

MATLAB信号处理工具中的各种IIR数字滤波器设计函数都是采用的双线性变换法。

题目1 设计一个Butterworth高通数字滤波器，通带边界频率为400Hz，阻带边界频率为300Hz，通带的波纹小于1dB，阻带衰减大于40dB,采样的频率为2000Hz。假设一个信号为x(t)=sin2*pi*f1*t+0.5cos2*pi*f2*t,其中f1=100Hz,f2=500Hz。试将原信号与通过该滤波器的输出信号进行比较。 解：其源程序的代码如下 Fs=2000; %采样频率

wp=400*2/Fs;ws=300*2/Fs; %根据采样频率将滤波器边界频率进行转换

Rp=1;Rs=40; %通带波纹与阻带衰减

Nn=128; %显示滤波器的频率特性的数据长度 [N,Wn]=buttord(wp,ws,Rp,Rs); %求得最小阶数和截止频率 [b,a]=butter(N,Wn,'high'); %设计Butterworth高通滤波器

3

figure(1)

[H,f]=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/db');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; n=0:127;

dt=1/Fs;t=n*dt; %时间序列 f1=100;f2=500; %输入信号频率

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 figure(2) subplot(2,1,1); plot(t,x); title('输入信号') y=filter(b,a,x); subplot(2,1,2); plot(t,y); title('输入信号') xlabel('时间/s') 其设计效果如下所示：

4

0振幅/db-100-200-3000100200300400500600频率/Hz70080090010005000相位/o-500-1000-15000100200300400500600频率/Hz7008009001000

图1 高通滤波器效果图

输入信号210-1-200.010.020.030.040.050.060.07输入信号0.50-0.500.010.020.030.04时间/s0.050.060.07

图2 信号分离效果图

5

2.3 冲激相应不变法的IIR数字滤波器

冲激响应不变法的设计原理是利用数字滤波器的单位抽样响应序列H(Z)来逼近模拟滤波器的冲激响应g(t)，按照冲激响应不变法原理，通过模拟滤波器的系统传递函数G(s)，可以直接求的数字滤波器的系统函数H(z)，其转换步骤如下：

（1）

利用=T，将p，s转换成p，，而p，s不变；

（2） （3）

求解低通模拟滤波器的传递函数G(s);

将模拟滤波器的传递函数G(s)转换为数字滤波器的传递函数H(z)。

题目二 设计模拟低通巴特沃斯滤波器，通带的波纹为Rp=1dB，通带上限角频率p=0.2，阻带下限角频率 s =0.3，阻带最小衰减s=15dB，根据该低通模拟滤波器，利用冲激响应不变法设计响应的数字低通滤波器，并且绘出设计后的数字滤波器的特性曲线。

解：其源程序的代码如下 wp=0.2*pi; ws=0.3*pi; Rp=1; As=15; T=1;

Rip=10^(-Rp/20); Atn=10^(-As/20); OmgP=wp*T; OmgS=ws*T;

[N,OmgC]=buttord(OmgP,OmgS,Rp,As,'s');%选取模拟滤波器的

6

阶数

[cs,ds]=butter(N,OmgC,'s');%设计所需要的模拟低通滤波器 [b,a]=impinvar(cs,ds,T);%应用脉冲响应不变法进行转换

[db,mag,pha,grd,w]=freqz_m(b,a);%求得相对、绝对频响及相位、群延迟

subplot(2,2,1); plot(w/pi,mag); title('幅频特性'); xlabel('w(pi)'); ylabel('|H(jw)|'); axis([0,1,0,1.1]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); set(gca,'YTickMode','manual','YTick',[0 Atn Rip 1]); grid

subplot(2,2,2); plot(w/pi,db);

title('幅频特性(db)'); xlabel('w(/pi)'); ylabel('dB');

axis([0,1,-40,5]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); set(gca,'YTickMode','manual','YTick',[-40 -As -Rp 0]); grid

subplot(2,2,3); plot(w/pi,pha/pi); title('相频特性'); xlabel('w(/pi)'); ylabel('pha(/pi)'); axis([0,1,-1,1]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); grid

subplot(2,2,4); plot(w/pi,grd); title('群延迟'); xlabel('w(/pi)');

7

ylabel('Sample'); axis([0,1,1,12]);

set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]); grid

%所使用的M文件函数：

function[db,mag,pha,grd,w]=freqz_m(b,a) [H,w]=freqz(b,a,500);%500点的复频响应 mag=abs(H);%绝对幅值响应

db=20*log(mag/max(mag));%相对幅值响应 pha=angle(H);相位响应

grd=grpdelay(b,a,w);%群延迟响应 其滤波图形如下所示：

幅频特性10.89130-1-15幅频特性(db)|H(jw)|0.1778000.20.30.5w(pi)相频特性1-4000.20.30.5w(/pi)群延迟110.50-0.5-100.20.30.5w(/pi)1dBSamplepha(/pi)1210864200.20.30.5w(/pi)1

图3 模拟低通巴特沃斯滤波器效果图

通过观察图，发现很好的满足了题目所给的要求。

8

3 FIR数字滤波器的设计

由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。但是这些特性是以牺牲线性相位为代价的，即用Butterworth、chelbchev和椭圆法设计的数字滤波器逼近理想的滤波器的幅度频率特性，得到的滤波器往往是非线性的。在许多的电子系统中，对幅度频率和线性相位特性都有较高的要求，所以IIR滤波器在这些系统中往往难以胜任。有限长的单位冲激响应（FIR）数字滤波器具有以下的优良特点：

（1） 可以在设计任意幅度频率特性滤波器的同时，保证精

度、严格的线性相位特性。

（2） FIR数字滤波器的单位冲激响应h(n)是有限长的，可

以用一个固定的系统来实现，因而FIR数字滤波器可以做成因果稳定的系统。 （3） 允许设计多通带系统。

3.1窗函数法设计FIR数字滤波器

窗函数法就是设计FIR数字滤波器的一种的方法。它在设计

FIR数字滤波器中很有重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular)，三角窗(Triangular window)，汉宁窗(Hanning window)，海明窗(Hamming window)，切比雪夫窗(Chebyshev window)，巴特利特窗(Bartlett window)。

9

3.2窗函数法设计步骤

根据给定的滤波器技术指标，选择滤波器长度N和窗函数n，使其具有最窄宽度的主瓣和最小的旁瓣。

窗函数法的设计步骤：

􀌈1） 给定理想频响函数 Hd(ejw)

􀌈2） 根据指标选择窗函数。确定窗函数类型的主要依据

过度带宽和阻带最小衰耗的指标；

􀌈3） 由Hd(ejw)求ha(n)，加窗得 h(n)=hd(n)(n) 􀌈4）检验，由h(n)求H(ejw)，求H(ejw)是否在误差容限之内。

题目3 根据下列技术指标，设计一个FIR数字低通滤波器： p=0.25*pi s =0.35*pi

Ap=0.25dB As= 50dB

选择一个适当的窗函数，确定单位冲激响应，绘出所涉

及滤波器的幅度响应。

首先根据窗函数的最小特性表，如下 窗函数 矩形窗 三角窗 汉宁窗 海明窗 第一旁瓣相对于主瓣衰主瓣宽 减/dB -13 4/N -15 -31 -41 8/N 8/N 8/N 12/N 阻带最小衰减/dB 21 25 44 53 74 布拉克曼窗 -57 10

凯塞窗 可调 可调 可调 可调 可调 切比雪夫窗 可调 其源程序代码如下： clear all; Wp=0.25*pi; Ws=0.35*pi;

tr_width=Ws-Wp; %过渡带宽度 N=ceil(6.6*pi/tr_width)+1; %滤波器长度

n=0:1:N-1; %理想低通滤波器的截止频率 Wc=(Ws+Wp)/2; %理想低通滤波器的单位冲激响应 hd=ideal_lp1(Wc,N); %海明窗

w_ham=(hamming(N))'; %截取得到实际的单位脉冲响应 h=hd.*w_ham; %计算实际滤波器的幅度响应 [db,mag,pha,w]=freqz_m2(h,[1]); delta_w=2*pi/1000;

Ap=-(min(db(1:1:Wp/delta_w+1))) %实际通带波纹 As=-round(max(db(Ws/delta_w+1:1:501))) %实际阻带波纹 subplot(221) stem(n,hd)

title('理想单位脉冲hd(n)') subplot(222) stem(n,w_ham) title('海明窗') subplot(223)

11

stem(n,h)

title('实际单位脉冲响应hd(n)') subplot(224) plot(w/pi,db) title('幅度响应(dB)') axis([0,1,-100,10]) grid

%所使用的M文件如下

function[db,mag,pha,w]=freqz_m2(b,a) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:1:501))'; w=(w(1:1:501))'; mag=abs(H);

db=20*log10((mag+eps)/max(mag));%相对幅值响应 pha=angle(H);%相位响应

运行结果： 34 Ap= 0.0316 As= 52

由此结果知，所设计的低通滤波器为II型滤波器，它的通带波纹和阻带波纹均满足设计要求。

其相应的图形如下所示：

12

理想单位脉冲hd(n)0.30.20.10-0.102040608000200.51海明窗406080实际单位脉冲响应hd(n)0.30.20.10-0.1020406080-1000-500幅度响应(dB)0.51

图4 窗函数滤波器效果图

3.3利用频率采样法设计FIR滤波器

窗函数设计FIR数字滤波器是从时域出发，把理想的滤波器的单位取样响应用合适的窗函数截短成为有限长度的h(n)，并使逼近理想的hd(n)，以实现所设计的滤波器的频率响应Hd(ejw)逼近理想滤波器的频率响应Hd(ejw)。

一个有限长的序列，如果满足频率采样定理，可以通过频率的有限个采样点的值被准确的得以恢复。

题目4 根据下列的技术指标，利用频率抽样方法1，设计一个1型的FIR数字高通滤波器：

Wp=0.8*pi Ws=0.7*pi Ap=1bB As=40bB

并绘出频率抽样响应以及所设计滤波器的单位冲激响应和幅度响

13

应图。

选择N=61, 则Wp位于k=24,Ws位于k=21，并令T1=0.1095,T2=0.598 解：程序设计如下 clear all; N=61;

T1=0.1095; T2=0.598;

alpha=(N-1)/2; l=0:N-1;

w1=(2*pi/N)*l;

Hrs=[zeros(1,22),T1,T2,ones(1,14),T2,T1,zeros(1,21)]; Hdr=[0,0,1,1];

wd1=[0,0.75,0.75,1]; k1=0:floor((N-1)/2); k2=floor((N-1)/2)+1:N-1;

angH=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)]; Hdk=Hrs.*exp(j*angH); h=real(ifft(Hdk,N));

[db,mag,pha,w]=freqz_m2(h,l) [Hr,ww,a,L]=hr_type1(h); subplot(221)

plot(w1/pi,Hrs,'.',wd1,Hdr) title('频率样本Hd(k):N=61') axis([0 1 -0.1 1.2]) subplot(222) stem(l,h)

title('实际单位脉冲响应h(n)') subplot(223)

plot(ww/pi,Hr,w1/pi,Hrs,'.') title('实际振幅响应H(w)') axis([0 1 -0.1 1.2]) subplot(224) plot(w/pi,db)

title('幅度响应(db)') axis([0 1 -80 10]

其运行的图形结果如下：

14

频率样本Hd(k):N=610.410.20-0.2000.5实际振幅响应H(w)10-200.5-40-60000.51-8001-0.40实际单位脉冲响应h(n)0.5204060幅度响应(db)0.51

图5 频率抽样法设计FIR滤波器

15