基于MATLAB的巴特沃斯滤波器

数字信号处理课程设计2015年 6 月25 日目录一．设计目的： (3)二．设计要求： (3)三．设计内容： (4)3.1选择巴特涡斯低通数据滤波器及双线性变换法的原因 (4)3.2巴特沃思低通滤波器的基本原理 (4)3.3双线性变换法原理 (5)3.4数字滤波器设计流程图 (7)3.5数字滤波器的设计步骤 (7)四．用matlab实现巴特沃斯低通数字滤波器的仿真并分析 (9)4.1巴特沃斯低通数字滤波器技术指标的设置 (9)4.2用matlab实现巴特沃斯低通数字滤波器的仿真 (9)4.3波形图分析： (12)五.总结与体会 (13)六．附录参考文献 (14)2一．设计目的：该课程设计是测控技术与仪器专业的必修课，开设课程设计的目的使学生掌握数字信号处理的基本概念和基本理论，能够利用辅助工具进行FIR和IIR数字滤波器的设计，进行一维信号的频谱分析，并进行仿真验证。

加强实践教学环节，加强学生独立分析、解决问题的能力，培养学生动手能力和解决实际问题的能力，实现宽口径教育。

(1)理解低通滤波器的过滤方法。

(2)进一步熟悉低通滤波器的基本应用。

(3)用仿真工具matlab软件对设计的滤波器进行软件和硬件仿真。

(6)将对仿真结果进行比较，从而检验滤波器滤波性能的准确性。

二．设计要求：地震发生时，除了会产生地震波，还会由地层岩石在断裂、碰撞过程中所发生的震动产生次声波。

它的频率大约在每秒十赫兹到二十赫兹之间（可以用11Hz和15Hz的两个信号的和进行仿真，幅度可以分别设定为1、2）。

大气对次声波的吸收系数很小，因此它可以传播的很远，而且穿透性很强。

通过监测次声波信号可以监测地震的发生、强度等信息，因为自然界中广泛存在着各种次声波，这就对地震产生的次声波产生了干扰（可以用白噪声模拟，方差为5），需要采取一定的处理方法，才能检测到该信号，要求设计检测方案；并处理方法给出具体的软件（可以以51系列单片机、STM32F407、TMS320F28335或TMS320F6745为例）。

假设地震次声波信号为x,输入x=sin(2*π*11*t)+2*sin(2*π*15*t)和伴有白噪声的合成信号，经过滤波器后滤除15Hz以上的分量，即只保留x=sin(2*π*11*t)+2*sin(2*π*15*t)的分量信号，来验证设计的滤波器是否达到了设计要求。

34 三．设计内容：3.1选择巴特涡斯低通数据滤波器及双线性变换法的原因（1）由于低通滤波器是组成其它滤波器的基础，故选用低通滤波器； （2）在当今社会，数字信号的应用越来越广泛，故选用数字信号；（3）巴特沃斯滤波器的特点是通频带的频率响应曲线最平滑并且应用范围最广，故选巴特沃斯型滤波器；（4）为了不使数字滤波器在ω=π附近产生频谱混叠，故选用双线性变换法。

3.2巴特沃思低通滤波器的基本原理巴特沃斯低通数字滤波器的幅度平方函数2a j H ）（Ω用下式表示Na c22)(11)j (H ΩΩ+=Ω 式中，N 称为滤波器的阶数。

当Ω=0时，1j H a =Ω）（；c Ω=Ω时，2/1j H a =Ω）（，c Ω是3dB 截止频率。

在c Ω=Ω附近，随Ω加大，幅度迅速下降。

幅度特性与Ω与N 的关系如图3.1所示。

幅度下降的速度与阶数N 有关，N 愈大，通带愈平坦，过渡带愈窄，过渡带与阻带幅度下降的速度愈快，总的频响特性与理想低通滤波器的误差愈小。

图3.1 巴特沃斯低通数字滤波器 图3.2 三阶巴特沃斯滤波器极点5 幅度特性与Ω与N 的关系 分布图Na a cj s s s 2)(11)(H )(H Ω+=- 复变量Ω+=j s σ，此式表示幅度平方函数有2N 个极点，极点ks 用下式表示：)21221(21()1(Nk j c Nk ec j s ++Ω=Ω-=π） （k =0,1,2,3….）2N 个极点等间隔分布在半径为c Ω的圆上（该圆称为巴特沃斯圆），间隔为N /πrad 。

例如N=3,极点间隔为π/3rad ，如图3.2所示。

为形成因果稳定的滤波器，2N 个极点中只取ｓ平面左半平面的的N 个极点构成Ha(s), 而右半平面的的N 个极点构成Ha(-s),Ha （s ）的表达式为)()(Ha 1k N k Ncs s s -Ω=∏-=为使设计公式和图表统一，将频率归一化。

巴特沃斯低通数字滤波器采用对3dB 截止频率c Ω归一化，归一化后的系统函数为∏=ΩΩ=Ω1-0k )c s -(1)(H N k cc s s a 令ληj p +=,c ΩΩ=/λ,λ称为归一化频率，p 称为归一化复变量，这样，巴特沃斯低通原型系统函数为∏==1-0)-(1Ga(p)N k kpp3.3双线性变换法原理双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。

为了克服多值映射的缺点，采用把整个s 平面频率压缩方法，6 将整个频率轴上的频率范围压缩到-π/T ～π/T 之间，再用sT e Z =转换到Z 平面上。

也就是说，第一步先将整个S 平面压缩映射到S1平面的-π/T ～π/T 一条横带里；第二步再通过标准变换关系sT e Z =将此横带变换到整个Z 平面上去。

这样就使S 平面与Z 平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。

映射关系如图3.3所示。

设Ha （s ）,Ω=j s ,经过非线性频率压缩后用^1)(s H a ，11Ω=j s 表示，这里用正切变换实现频率压缩： )21tan(21T T Ω=Ω图3.3 双线性变换的映射关系式中，T 为采样间隔，当1Ω从-π/T 经过0变化到π/T 时，Ω由-∞经过0变化到+∞，实现了s 平面上整个虚轴完全压缩到1s 平面上虚轴的+π/T 之间的转换。

即T j Tj T j T j tT j T j e e T ee e e T j 11111111222/2/2/2/Ω-Ω-ΩΩΩ-Ω+-=+-=Ω代入Ω=j s ,11Ω=j s ,得到Ts Ts e e T s 11112--+-=再通过T s e z1=从1s 平面转换到z 平面，得到7 11112--+-=z z T s s Ts T z -+=22上式是S 平面与Z 平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换。

双线性变换法与冲激响应不变法相比，其主要的优点是避免了频率响应的混叠现象，虽然在线性方面有些欠缺，但是可以通过频率的预畸来加以校正且计算比冲激响应不变法方便，实现起来比较容易，所以，本设计选择用双线性变换法设计巴特沃斯低通滤波器。

3.4数字滤波器设计流程图3.5数字滤波器的设计步骤数字滤波器的设计步骤：根据数字滤波器的技术指标先设计过渡模拟滤波器得到系统函数Ha(s),然后将Ha(s)按某种方法（本实验采用双线性变换法）转换成数字滤波器的系统函数H （z ）。

具体为：（1）确定巴特沃斯数字低通滤波器的技术指标：通带边界频率ωp,阻带截止频率ωs,通带最大衰减аp,阻带最小衰减аs 。

（2）将数字滤波器的技术指标转换为模拟滤波器的技术指标。

这里指ωp 和ωs 的变换而аp 和аs 保持不变。

本题采用双线性变换法，其转换公式为：2tan 2p T p ω=Ω8 2tan 2s T s ω=Ω （3）根据技术指标Ωp 、Ωs 、ωp 和ωs 用下面公式求出滤波器的阶数。

p sΩΩ=sp λ 11011010/10/--=p s sp k αα sp g sp g l k l λ=N（4）根据N 由表3.1求出归一化极点k p 和归一化低通原型系统函数Ga(p)。

表3.1巴特沃斯归一化低通滤波器参数（5）将Ga （p ）去归一化，将csp Ω=代入Ga （p ），得到实际的滤波器系统函数： cs p p G s Ha Ω==)()(这里Ωc 为3dB 截止频率。

9 （6）用双线性变换法将模拟滤波器Ha(s)转换成数字低通滤波器系统函数H(z)。

转换公式为：)()(s Ha z H =s=11112--+-zz T 四．用matlab 实现巴特沃斯低通数字滤波器的仿真并分析4.1巴特沃斯低通数字滤波器技术指标的设置 通带截至频率ωp=15Hz, αp =1dB 阻带截至频率ωs=20Hz αs =30dB 采样频率为fs=1000Hz4.2用matlab 实现巴特沃斯低通数字滤波器的仿真 Matlab 程序如下： clear all;%模拟地震信号，频率是11hz 和15hz fs=1000;dt=1/fs; f1=11;f2=15; n=500;t=[0:n-1]*dt; %时间序列 x=sin(2*pi*f1*t)+2*sin(2*pi*f2*t); %信号 figure(1); subplot(511);plot(t,x); %显示原始信号 title('模拟地震信号'); %白噪声信号%rand 函数用来产生均值0.5，方差约为1/12，幅值在0~1的伪随机数 %修改为均值为0，方差为5的白信号。

p=5; %u1=rand(1,n);u1_mean=mean(u1);u1_var=var(u1);u=u1-u1_mean;u=u*sqrt(p/u1_var); %白噪声信号subplot(512);plot(u(1:100));title('均匀分布白噪声');%%%%地震信号和白噪声叠加y1=x+u; %叠加白噪声subplot(513);plot(t,y1);title('地震信号和白噪声叠加');%FIR带通%m=20;%f=[0 0.001 0.0015 0.004 0.005 1];%a=[0 0 1 1 0 0];%BB=firls(m,f,a);%bb=fftfilt(BB,y1);%subplot(313);%plot(t,bb);%fir低通%m=60;%f=[0 0.03 0.04 1];%a=[1 1 0 0];%BB=firls(m,f,a);10%bb=fftfilt(BB,y1);%subplot(313);%plot(t,bb);%B=fir1(45,0.025);%bb=fftfilt(B,y1);%subplot(313);%plot(t,bb);%巴特沃斯wp=2*20/fs;ws=2*30/fs;Rp=1;As=30;[N,wc]=buttord(wp,ws,Rp,As);[B,A]=butter(N,wc);bb=filter(B,A,y1);[H,W]=freqz(B,A);subplot(514);plot(W,abs(H));title('巴特沃斯幅频特性');subplot(515);plot(t,bb);title('巴特沃斯滤波');用matlab滤波前后的信号波形变化如图4.1所示：111210输入信号00.51 1.52 2.53 3.5012低通滤波器00.20.40.60.81230Hz图4.1 用matlab 滤波前后的信号波形变化4.3波形图分析：由技术指标得：设计的巴特沃斯低通数字滤波器为15Hz 以内的信号能通过，而高于15Hz 的信号将通不过滤波器。