基于MATLAB的IIR数字带通滤波器的设计.doc

1、基于MATLAB的IIR数字带通滤波器的设计设计一:利用冲击响应不变法设计一个切比雪夫带通滤波器,设计的指标为: 通带边缘频率 w1 = 0.4pi,w2 = 0.8pi 阻带边缘频率 w3 = 0.3pi,w4 = 1.0pi 通带波纹 = 0.5dB阻带衰减 = 20 dB解: 1.基本原理 滤波器的传递函数的一般形式为:n 当M=N, H(z): N阶IIR系统+(M-N)阶的FIR系统,n 以上两种表示等价，部分分式形式和零极点增益形式n IIR系统的逼近, 就是找到滤波器的系数ak, bk，或者是系统的零极点和增益(z,p,k)。切比雪夫I型的幅度平方函数为: 的特点如下:(1) 当

2、=0时,N为偶数时,=,当N为奇数时，=0.(2) 当=c时， 即此时所有的幅度函数曲线都经过点, c即为切比雪夫滤波器的通带截止频率.(3) 在通带内,在1之间等波纹地起伏。(4) 在通带外,随着的增大, 迅速单调的趋近于零。该滤波器在通带内具有等波纹起伏特性,在阻带内则单调下降且具有更大的衰减.相比于巴特沃斯滤波器,阶数N较小. 2 冲击响应不变法 冲击响应不变法是使数字滤波器的单位冲击响应序列h(n)模仿模拟滤波器的单位冲击响应ha(t)。将模拟滤波器的单位冲击响应加以等间隔抽样,使h(n)正好等于ha(t)的抽样值，即满足:h(n)= ha(nT)既有: 冲击响应不变法是将模拟滤波器的

3、s平面变换成数字滤波器的z平面。冲击响应不变法使得数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，也就是时域逼近良好。一个线性相位的模拟滤波器通过脉冲响应不变法得到的仍然是一个线性相位的数字滤波器。 但是，由于任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生频率响应的混叠失真。这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。3设计流程: (1)确定数字滤波器的技术指标。 (2)将数字滤波器的技术指标转换成相

4、应的模拟滤波器的技术指标。 (3)按照模拟低通滤波器的技术指标设过渡模低通滤波器。 (4)用脉冲响应不变法，将模拟滤波转换成数字滤波器 4.matlab实现主要步骤:(1)先由数字频率转化为模拟的频率,由于采用脉冲响应不变法,所以: Omega1 =fs*w1; Omega2 = fs*w2; Omega3 =fs*w3; Omega4 = fs*w4;采样频率fs = 2000HZ 由函数 N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s); z,p,k = cheb1ap(N,Rp); 求出模拟低通滤波器。(2)再通过bt,at = lp2bp(bb,aa,Wo,Bw)函数实现由模拟

5、低通滤波器向带通滤波器的转换. 其中:Wo = sqrt(Omega1*Omega2); Bw = Omega2 - Omega1; 进行频率变换后实现了带通模拟滤波器(3): 采用脉冲响应不变法:利用函数: bz, az=impinvar(bt,at,fs) 将上述带通模拟滤波器转换为数字带通滤波器. 其中fs为采样频率(4):验证： 预期输出信号: uzs=sin(2*pi*451*nn) 加入噪声后的信号为: uz = 0.5*cos(2*pi*21*nn) + uzs+ 2*sin(2*pi*1000*nn)观察经过该滤波器后的输出波形。 其中Fn = 100; dn= 1/fs;n=

6、0:Fn-1; nn= n*dn; 5.matlab实现结果: (1)根据滤波器的设计指标,首先实现一个模拟的低通滤波器,其曲线如图一所示: 图一：模拟低通滤波器 (2) 将设计的模拟低通滤波器进行频率变换,变换成模拟的带通滤波器,其幅频,相频曲线如图2所示: 图二:模拟带通滤波器 (3)冲击响应不变法实现数字带通滤波器的带通曲线如图三所示: 其幅频相频曲线如图四所示: 图四:带通滤波器的幅频和相频曲线(5) 为了对设计的数字带通滤波器的性能进行验证,本设计加入了噪声信号,如图五所示即为原来理想的信号波形 uzs=sin(2*pi*451*nn)(6) 如图六所示即为加入噪声信号后的波形曲线。

7、 噪声信号的曲线方程为: uz=0.5*cos(2*pi*21*nn)+uzs+2*sin(2*pi*1000*nn)(7) 如图七所示即为通过滤波器后输出的波形 图五:理想波形图六:噪声信号图七:滤波后输出波形Matlab程序代码:w1 = 0.4*pi,w2 = 0.8*pi; %digital signal freqRp=0.5;Rs = 20; w3 = 0.3*pi,w4=1*pi; fs = 2000;Omega1 =fs*w1; Omega2 = fs*w2;Omega3 =fs*w3;Omega4 = fs*w4;f1 =Omega1/2/pi;f2 = Omega2/2/pi

8、;f3 =Omega3/2/pi;f4 = Omega4/2/pi;Wp = Omega1,Omega2;Ws = Omega3,Omega4;N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s);z,p,k = cheb1ap(N,Rp);bs,as = zp2tf(z,p,k);n = 0:0.01:2;hs,ws=freqs(bs,as,n);figure;plot(ws,abs(hs).2); %low pass filtergrid;bb=k*real(poly(z);aa=real(poly(p);Wo = sqrt(Omega1*Omega2);Bw = Omega2 - O

9、mega1;bt,at = lp2bp(bb,aa,Wo,Bw);h,w = freqs(bt,at);figure;subplot(2,1,1);plot(w/2/pi,(abs(h).2);axis(100 1500 0 1);grid;subplot(2,1,2);plot(w/2/pi,angle(h)*180/pi);grid;H = tf(bt,at);bz, az=impinvar(bt,at,fs);hz,wz=freqz(bz,az,fs);figure;plot(wz/pi,abs(hz);grid;figure;freqz(bz,az);Fn = 100;dn= 1/fs

10、;n=0:Fn-1;nn= n*dn;uzs=sin(2*pi*451*nn);uz = 0.5*cos(2*pi*21*nn) + uzs+ 2*sin(2*pi*1000*nn);figure;plot(nn,uz);grid;axis(0 0.05 0 2);ylabel(uz);Y = filter(bz,az,uz);figure;plot(nn,uzs);axis(0 0.05 0 1.5);xlabel(预输出信号波形);grid;figure;plot(nn,Y);axis(0 0.05 0 1.5);xlabel(实际输出波形);grid; 设计二: 设计一Kalman滤波器

11、的设计1.概要:Kalman滤波理论是Wiener滤波理论的发展，具有如下特点(1):数学公式用状态空间概念描述；(2):它的解释递推计算的，与Wiener滤波器不同，Kalman滤波器是一种自适应滤波器。对于一般的时变系统,系统状态方程和观测方程可表示为一下形式: (1) 状态方程(state equation) 状态向量： 状态转移矩阵： 状态噪声输入矩阵： 系统状态噪声： 其中噪声 为均值为零的高斯白噪声(2) 观测方程(measurement equation) 观测向量: 观测矩阵: 观测噪声: 其中噪声 为均值为零的高斯白噪声 2.应用: 假设一理想质点的运动轨迹是一理想的正弦波.

12、如:y=5*sin (pi*t*0.5) 但是当质点受到外部环境的干扰后就会产生一定的偏移,利用函数 noise=randn(1,200) 产生随机的信号进行模拟,所以实际的质点偏离中心位移的曲线为y_true= y + noise.可令质点的运动模型为y = sin (x),质点做匀速直线运动。令其经过的位移为s(n),速度为v(n),所以有: 对于质点的干扰可以视为质点产生的加速度,从而有: 其中:a表示质点受到的加速度 T表示采样的间隔从而状态方程可以表示为:令x(n)=(s(n),v(n)T, ,v1(n-1)=a(n-1) 状态方程化为: 其中状态转移矩阵F(n,n-1)=观测方程z

13、(n)=C(n)x(n) + V2(n)其中C(n)=(1 0); V2(n)=noise3.matlab实现: (1)理想质点的运动轨迹如图一所示: 图一:理想质点运动轨迹(2)实际的质点的运动轨迹,在matlab中用加上均值为零的高斯白噪声进行模拟，波形如图二所示: 图二：实际的质点运动轨迹 (3) 经过Kalman滤波器后的输出的波形如图三所示: 图三: 输出质点的位移4.matlab程序代码%-产生正弦波信号-%T=0.1;t=0:T:20-0.1;y=5*sin(pi*t*0.5);figure;plot(t,y);grid;axis(0,20,-8,8);xlabel(时间/s);

14、ylabel(质点位移/mm);title(理想质点运动轨迹);%-添加噪声-%noise=randn(1,200);noise = randn(1,200);noise=noise-mean(noise);noise=sqrt(var(noise)*noise/sqrt(var(noise);Z=y+noise;figure;plot(t,Z);grid;axis(0,20,-8,8);xlabel(时间/s );ylabel(质点位移/mm);title(实际的质点运动轨迹);%-建立系统模型-%X=zeros(2,200);X(:,1)=0,1;F=1,T;0,1; V=1/2*(T)2

15、 T;%V=T T; how to get the value of matix V%V=1 1;C=1,0;P=0,0;0,1;%Q1=(0.25)2; %Q2=(0.25)2; Q1=cov(noise); Q2=cov(noise); %-卡尔曼算法-%for n=1:200 Kg=P*C/(C*P*C+Q2); X(:,n)=X(:,n)+Kg*(Z(:,n)-C*X(:,n); X(:,n+1)=F*X(:,n); P=(eye(2,2)-Kg*C)*P; P=F*P*F+V*Q1*V;end%-%figure t=0:0.1:20;Y=X(1,:);plot(t,Y);grid;axis(0,20,-8,8);xlabel(时间/s );ylabel(质点位移/mm );title(卡尔曼滤波后输出质点最优化运动位移); .