clear;close all;clc;
fs=100;
delta_t=1/fs;n=0:1:255;
fp1=10;
fp2=30;
w1=2*pi*fp1/fs;
w2=2*pi*fp2/fs;
theta1=(w1+w2)/2;
theta2=(w2-w1)/2;
k=cot(theta2)*tan(2*pi*10/fs);
gama=cos(theta1)/cos(theta2);
beta=(k-1)/(k+1);
alpha=(2*gama*k)/(k+1);
c=2*pi*10/fs;
d1=c^2+sqrt(2)*c+1;
d2=c^2-sqrt(2)*c+1;
b0=c^2/d1;
b1=2*b0;
b2=b0;
a1=2*(c^2-1)/d1;
a2=d2/d1;
bp0=b0-b1*beta+b2^2*beta^2;
bp1=-2*alpha*b0+alpha*b1+alpha*beta*b1-2*alpha*beta*b2;
bp2=alpha^2*b0+2*b0*beta-b1-alpha^2*b1-b1*beta^2+alpha^2*b2+2*beta*b2;
bp3=-2*alpha*beta*b0+alpha*b1+alpha*beta*b1-2*alpha*b2;
bp4=beta^2*b0-beta*b1+b2;
ap0=1-a1*beta+a2*beta^2;
ap1=-2*alpha+alpha*a1+alpha*beta*a1-2*alpha*beta*a2;
ap2=alpha^2+2*beta-a1-alpha^2*a1-beta^2*a1+alpha^2*a2+2*beta*a2;
ap3=-2*alpha*beta+alpha*a1+alpha*beta*a1-2*alpha*a2;
ap4=beta^2-a1*beta+a2;
B=[bp0,bp1,bp2,bp3,bp4];
A=[ap0,ap1,ap2,ap3,ap4];
AA=round(256.*A);
BB=round(256.*B);
x2=sin(2*pi*5*n*delta_t)+sin(2*pi*20*n*delta_t)+sin(2*pi*35*n*delta_t);
y2=filter(B,A,x2);
figure(2)
freqz(B,A);
title('滤波器幅频相应');
figure(1)
subplot(3,1,1);
plot(n,x2);
grid on;
title('原始信号')
subplot(3,1,2);
plot(n,y2);
grid on;
title('滤波后信号')
subplot(3,1,3);
plot(n,x2-y2);
title('滤掉信号')
grid on;
