batewosidaitong.m

基于MATLAB巴特沃思滤波器的研究

clear%清除内存
f1=300;f2=500;fs1=100;fs2=600;%实际频率
fs=2000;%采样频率
t=1/fs;%采样时钟
wp1=2*f1*pi/fs;%数字滤波器的通带带下限角频率
wp2=2*pi*f2/fs;%数字滤波器的通带上限角频率
ap=3;   %带边频率的衰减
ws1=2*pi*fs1/fs;%数字滤波器的下阻带角频率
ws2=2*pi*fs2/fs;%数字滤波器的上阻带角频率
as=18;% 阻带最大的摔减
wc1=(2/t)*tan(wp1/2);wc2=(2/t)*tan(wp2/2);%把数字滤波器的技术指标转模滤波器技术指标频率预畸变
wr1=(2/t)*tan(ws1/2);wr2=(2/t)*tan(ws2/2);%...............................................
wo=sqrt(wc1*wc2);%求中心频率
    BW=wc2-wc1;%求带宽
nwr1=(((wr1^2)-(wo^2))/(BW*wr1));%以带宽作为参考做归一化并把带通滤波器归一化转为低通归一化
nwr2=(((wr2^2)-(wo^2))/(BW*wr2));...........................................................
nwc1=(((wc1^2)-(wo^2))/(BW*wc1));.......................................................
nwc2=(((wc2^2)-(wo^2))/(BW*wc2)); ......................................................
if(abs(nwr1)-abs(nwr2))>0   %判断
    nowr=abs(nwr2);
else           %否则
    nowr=abs(nwr1);
end             %取绝对值最小的
nowrc=1;   %低通归一化截至频率为1
n2=(log10(sqrt(10.^(abs(as)*0.1)-1)./(10.^(abs(ap)*0.1)-1)))./(log10(nowr));%计算模拟滤波器的阶次
n=ceil(n2);%取想正无穷取整
[z,p,k]=buttap(n);%设计模拟滤波器原型零点极点由于没有零点因此为空阵
[b,a]=zp2tf(z,p,k);%求模拟滤波器原型的系数
[b1,a1]=lp2bp(b,a,wo, BW);%将低通模拟滤波器原型转换为中心频率带宽的带通模拟滤波器的系数
[b2,a2]=bilinear(b1,a1,fs);%将模拟带通滤波器的系数转为数字滤波器的系数双线性法
w=[0:512]*pi/512;
hz=freqz(b2,a2,w);%计算向量w指定数字频率点的频率响应画出数字滤波器 的幅频特性
figure(1);%第1仿真图
subplot(3,1,1)%分3行第1图
plot(w/pi,20*log10(abs(hz)));grid;%归一化后绘图
xlabel('w/pi');
ylabel('分贝');
title('数字带通滤波器幅频响应');
line([0.3,0.3],[-25,0]);
line([0.5,0.5]  ,[-25,0]);
axis([0,1,-25,5])%坐标的范围
ph2=unwrap(angle(hz));%相频曲线
subplot(3,1,2);
plot(ph2);grid on;
xlabel('fs');
ylabel('角度');
title('数字带通相频响应');
w5=[0:800*2*pi];
h4=freqs(b1,a1,w5);%画出模拟滤波器幅频特性
subplot(3,1,3)
plot(w5/2/pi,abs(h4));grid on;
title('模拟滤波器');
figure(2);%第2仿真图
n1=0:500;%取点数
x1=cos(0.1*pi*n1)+sin(0.8*pi*n1);
w1=randn(size(x1));%产生正态为零白噪声
x=x1+w1;
subplot(2,1,1)
stem(x,'.');
xlabel('n');
ylabel('x(n)');
title('输入信号');
axis([0,60,-4,4])
y=filter(b2,a2,x);%对x向量的数据进行直接滤波即是差分方程求解
subplot(2,1,2)
stem(y,'.');
xlabel('n1');
ylabel('y(n)');
title('滤波后信号');
axis([0,60,-4,4])