数字信号处理结课论文.docx
《数字信号处理结课论文.docx》由会员分享,可在线阅读,更多相关《数字信号处理结课论文.docx(8页珍藏版)》请在冰豆网上搜索。
数字信号处理结课论文
基于MATLAB的数字滤波器设计
摘要
数字滤波器的实现是数字信号处理中的重要组成部分,设计过程较为复杂,牵涉到模型逼近、指标选择、计算机仿真、性能分析及可行性分析等一系列的工作,本文从设计原理以及数学软件matlab出发阐述数字滤波器的设计原理与方法。
关键词:
MATLAB,数字信号处理,数字滤波器
1绪论
数字滤波器是数字信号处理的重要应用,21世纪数字滤波器及其相关技术广泛的应用于通信、电子、自动控制等领域,是一种有效的抑制噪音、提取有用信号的方法。
它本身可以用硬件实现,也可以通过软件来实现,还可以通过专用的DSP处理器配合相应的软件,即软硬结合的方式来实现。
数字滤波器可以分为有限冲击响应(FIR)和无限冲激响应(IIR)两种。
通过MATLAB程序,实现输入相应技术指标及滤波器模型,输出相应数字滤波器的参数的功能。
2无限长单位脉冲响应滤波器(IIR)的设计
根据模拟滤波器设计数字滤波器,就是通过已知的模拟滤波器系统的系统函数H(s)来设计数字滤波器的系统函数H(z),主要是通过脉冲响应不变法,或双线性变换法完成S平面到Z平面的转换。
通过典型的模拟滤波器(诸如:
巴特沃斯滤波器、切比雪夫滤波器等)可以实现一定参数要求的数字滤波器。
根据已有的数字滤波器设计不同参数或者不同频带通断类型的数字滤波器。
例如已知数字低通滤波器的模型,通过变量代换得到不同截止频率的数字低通滤波器,或通过已知低通滤波器的模型设计高通、高阻、带通、带阻滤波器,这里主要是通过来完成相应的变量代换来实现滤波器类型的变换和参数的变换。
3有限长单位脉冲响应滤波器(FIR)的设计
IIR滤波器可用于较少的阶数达到所要求的幅度特性,且实现时所需的运算次数及存储单元都很少,十分适合于对于相位特性没有严格要求的场合,如果对相位特性有要求,这时选用FIR滤波器较好。
3.1窗函数法
从时域出发,把理想的窗口函数hd(n)截取成有限长的,以此h(n)来逼近理想的窗口函数hd(n),从而频率响应H(jw)也近似于理想的频率响应Hd(jw)。
主要窗函数有:
矩形窗,汉宁(Hanning)窗,汉明(Hamming)窗,布莱克曼(Blackman)窗,凯塞(Kaiser)窗。
3.2频率采样法
从频率出发,对理想的频率响应Hd(ejw)加以等间隔采样,以此Hd(k)实际FIR滤波器的频率特性的离散样本H(k)。
4利用MATLAB具体实现数字滤波器
4.1MATLAB中IIR数字滤波器设计的设计函数
信号处理工具箱提供的IIR经典设计方法是基于经典的低通滤波器到具有相同性能指标的数字滤波器的变换。
基本原理就是先根据滤波器的技术指标设计出相应的模拟滤波器,然后再将设计好的模拟滤波器变换为满足指标的数字滤波器。
利用表1的设计函数,可以很容易地产生任何阶数的高通、低通、带通、带阻数字滤波器。
表1滤波器设计函数
滤波器类型
设计函数
Bessel(用于模拟滤波器)
[b,a]=besself(n,Wn,options)
[z,p,k]=besself(n,Wn,options)
[A,B,C,D]=besself(n,Wn,options)
Butterworth
[b,a]=butter(n,Wn,options)
[z,p,k]=butter(n,Wn,options)
[A,B,C,D]=butter(n,Wn,options)
ChebyshevTypeI
[b,a]=cheby1(n,Wn,options)
[z,p,k]=cheby1(n,Wn,options)
[A,B,C,D]=cheby1(n,Wn,options)
ChebyshevTypeII
[b,a]=cheby2(n,Wn,options)
[z,p,k]=cheby2(n,Wn,options)
[A,B,C,D]=cheby2(n,Wn,options)
椭圆滤波器
[b,a]=ellip(n,Wn,options)
[z,p,k]=ellip(n,Wn,options)
[A,B,C,D]=ellip(n,Wn,options)
表中,n为滤波器阶次,Wn为滤波器的归一化截止频率;函数默认为低通或带通滤波器。
b,a分别为滤波器传递函数的分子和分母的系数向量;z,p,k分别为滤波器的零点、极点和增益。
options为滤波器类型参数,high为高通滤波器,截止频率为Wn;stop为带阻滤波器,截止频率为Wn=[W1,W2]。
4.2脉冲响应不变法设计数字滤波器
调用格式:
[bz,az]=impinvar(b,a,Fs),再给定模拟滤波器参数b,a和取样频率Fs的前提下,计算数字滤波器的参数。
两者的冲激响应不变,即模拟滤波器的冲激响应按Fs取样后等同于数字滤波器的冲激响应。
4.3利用双线性变换法设计数字ButterWorth滤波器
调用格式:
[bz,az]=bilinear[b,a,Fs],根据给定的分子b,分母系数a和取样频率Fs,根据双线性变换将模拟滤波器变换成离散滤波器,具有分子系数向量bz和分母向量az。
4.4IIR设计实例
利用脉冲响应不变法设计一个巴特沃斯数字低通滤波器,要求:
ωp=0.25π,Rp=1dB,ωs=0.4π,As=15dB,滤波器采样频率Fs=2000Hz。
程序如下:
wp=0.25*pi;
ws=0.4*pi;
Rp=1;As=15;
ripple=10^(-Rp/20);%滤波器的通带衰减对应的幅度值
Attn=10^(-As/20);%滤波器的阻带衰减对应的幅度值
%转换为模拟滤波器的技术指标
Fs=2000;T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs;
%模拟原型滤波器的计算
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,’s’)%计算阶数n和截止频率
[z0,p0,k0]=buttap(n);%设计归一化的巴特沃斯模拟滤波器原型
ba1=k0*real(poly(z0));%求原型滤波器的系数b
aa1=real(poly(p0));%求原型滤波器的系数a
[ba,aa]=lp2lp(ba1,aa1,Omgc);%变换为模拟低通滤波器
%用脉冲响应不变法计算数字滤波器的系数
[bd,ad]=impinvar(ba,aa,Fs)
[C,B,A]=dir2par(bd,ad)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H));
subplot(2,2,1);plot(w/pi,abs(H));
ylable(‘|H|’);title(‘幅度响应’);axis([0,1,0,1.1]);
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);
set(gca,’YTickMode’,’manual’,’YTick’,[0,Attn,ripple,1]);grid
subplot(2,2,2);plot(w/pi,dbH);title(‘幅度响应(dB)’);
ylabel(‘dB’);xlable(‘频率’);axis([0,1,-40,5]);
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);
set(gca,’YTickMode’,’manual’,’YTick’,[-50,-15,-1,0]);grid
subplot(2,2,4);zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title(‘零极点图’);
自定义dir2par函数:
functionI=cplxcomp(p1,p2)
I=[];
forj=1:
length(p2)
fori=1:
length(p1)
if(abs(p1(i)-p2(j))<0.0001)
I=[I,i];
end
end
end
I=I';
function[C,B,A]=dir2par(num,den)
M=length(num);
N=length(den);
[r1,p1,C]=residuez(num,den);
p=cplxpair(p1,10000000*eps);
I=cplxcomp(p1,p);
r=r1(I);
K=floor(N/2);
B=zeros(K,2);
A=zeros(K,3);
ifK*2==N;
fori=1:
2:
N-2;
Brow=r(i:
1:
i+1,:
);
Arow=p(i:
1:
i+1,:
);
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix(i+1)/2,:
)=real(Brow);
A(fix(i+1)/2,:
)=real(Arow);
end
[Brow,Arow]=residuez(r(N-1),p(N-1),[]);
B(K,:
)=[real(Brow),0];
A(K,:
)=[real(Arow),0];
else
fori=1:
2:
N-1;
Brow=r(i:
1:
i+1,:
);
Arow=p(i:
1:
i+1,:
);
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix(i+1)/2,:
)=real(Brow);
A(fix(i+1)/2,:
)=real(Arow);
end
end
程序运行结果如下:
n=5
Omgc=103.2016
bd=0.00720.03620.07250.07250.03620.0072
ad=1.0000-1.94341.9680-1.07020.3166-0.0392
sos=
1.00000.995601.0000-0.31930
1.00002.00721.00721.0000-0.69840.2053
1.00001.99720.99731.0000-0.92570.5976
g=0.0072
频率特性如图1所示:
图1频率特性曲线
4.5MATLAB中FIR数字滤波器的设计的函数
在MATLAB中产生窗函数十分简单:
根据长度n产生一个相应类型窗函数的类型:
(1)矩形窗(RectangleWindow):
调用格式:
w=boxcar(n)。
(2)三角窗(TriangularWindow):
调用格式:
w=triang(n)。
(3)汉宁窗(HanningWindow):
调用格式:
w=hanning(n)。
(4)哈明窗(HammingWindow):
调用格式:
w=hamming(n)。
(5)布莱克曼窗(BlackmanWindow):
调用格式:
w=blackman(n)。
(6)凯塞窗(KaiserWindow):
调用格式:
w=kaiser(n,beta)。
基于窗函数的FIR滤波器设计利用MATLAB提供的函数firl来实现。
调用格式:
firl(n,Wn,’ftype’,Window),n为阶数,Wn是截止频率(如果输入是形如[W1W2]的矢量时,本函数将设计带通滤波器,其通带为W1<ω格式中:
n为FIR滤波器的阶数,对于高通、带阻滤波器,n取偶数。
Wn为滤波器的截止频率,0~1;对于带通、带阻滤波器,Wn=[W1,W2],且W1‘ftype’为滤波器类型,默认为低通滤波器或带通滤波器。
其选择如下:
‘high’为高通滤波器,’stop’为带阻滤波器,’DC-1’使多带的第一频带为通带,’DC-0’使多带的第一频带为阻带。
Window为窗函数,列向量,其长度为n+1。
4.6FIR设计实例
用矩形窗设计一个FIR数字低通滤波器,要求:
N=64,截止频率wc=0.4,描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。
程序如下:
wc=0.4*pi;
N=64;n=0:
N-1;
hd=ideal_lp(wc,N);
windows=(boxcar(N))';
b=hd.*windows;
[H,w]=freqz(b,1);
dbH=20*log10(abs(H)+eps)/max(abs(H));
subplot(2,2,1),stem(n,hd);
axis([0,N,1.1*min(hd),1.1*max(hd)]);title('理想脉冲响应');
xlabel('n');ylabel('hd(n)');
subplot(2,2,2);stem(n,windows);
axis([0,N,0,1.1]);title('窗函数特性');
xlabel('n');ylabel('wd(n)');
subplot(2,2,3);stem(n,b);
axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');
xlabel('n');ylabel('h(n)');
subplot(2,2,4);plot(w/pi,dbH);
axis([0,1,-80,10]);title('幅度频率响应');
xlabel('频率(单位:
\pi)');ylabel('H(e^{j\omega})');
set(gca,'XTickMode','manual','XTick',[0,wc/pi,1]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
程序运行结果如图2所示:
图2FIR数字低通滤波器
总结
本文结合数字信号处理技术与计算机技术,从理论上分析了各种数字滤波器的设计思路,以及介绍了利用MATLAB进行各种数字滤波器的设计的基本方法。
利用MATLAB提供的各种库函数,及仿真工具大大的简化了设计步骤,使设计更方便、灵活。
体现了电脑辅助设计软件的优越性,为以后的数字滤波器设计学习与工作提供了新思路。
参考文献
[1]chi-Tsongchen.数字信号处理频谱计算与滤波器设计[M].北京:
电子工业出版社,2002.
[2]方勇.数字信号处理——原理与实践[M].北京:
清华大学出版社,2005.
[3]吴镇杨.数字信号处理[M].北京:
高等教育出版社,2005.
[4]刘波,文忠,曾涯.MATLAB信号处理[M].北京:
电子工业出版社,2002.