用双线性变换法设计原型低通为巴特沃兹型的数字IIR低通滤波器

clc;
clear all
Rp = 1; % 通带衰减/dB
Rs = 40; %阻带衰减/dB
OmegaP1_1=400; % 通带边界频率
OmegaS1_1=500; % 阻带边界频率
Fp=2000; % 抽样频率
Wp1=2*pi*OmegaP1_1/Fp; % 通带频率模数转换
Ws1=2*pi*OmegaS1_1/Fp; % 阻带频率模数转换
OmegaP1=2*Fp*tan(Wp1/2); % 通带频率非线性变换
OmegaS1=2*Fp*tan(Ws1/2); % 阻带频率非线性变换
Eta_P=OmegaP1/OmegaP1; %归一化
Eta_S=OmegaS1/OmegaP1;%归一化
Lemta_P_EquivalentLowPass=Eta_P;
Lemta_S_EquivalentLowPass=Eta_S; % 估计滤波器阶数
[N,Wn]=buttord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass, Rp, Rs,'s'); %设计滤波器
[num1,den1]=butter(N,Wn,'s');
[num2,den2]=lp2lp(num1,den1,OmegaP1);
[num,den]=bilinear(num2,den2,Fp);
[z,p,k]=tf2zp(num,den) %显示传输函数
w = 0:pi/255:pi;
h = freqz(num,den,w);
g = 20*log10(abs(h)); % 绘制增益响应
figure;
plot(w/pi,g);grid %绘制巴特沃兹低通滤波器幅频特性
axis([0 1 -60 5]);
xlabel('\omega /\pi');
ylabel('增益/ dB');
title('巴特沃兹低通滤波器');
figure;
zplane(z,p);%绘制传输函数零极点
title('传输函数的零极点');
M=100;
n=0:M;
x1=sin(2*pi*0.15*n);
x2=sin(2*pi*0.3*n);
x=x1+x2;
figure;
subplot(2,2,1)%绘制x1的波形
plot(x1);
grid on;
axis([0,100,-2,2]);
xlabel('t');
ylabel('x1(t)');
title('x1的波形');
subplot(2,2,2)%绘制x1的波形
plot(x2);
grid on;
axis([0,100,-2,2]);
xlabel('t');
ylabel('x2(t)');
title('x2的波形');
subplot(2,2,3)%绘制输入x的波形plot(x);
grid on;
axis([0,100,-2,2]);
xlabel('t');
ylabel('x(t)');
title('输入信号x的波形') hanshu=impz(num,den,600);
y=juanji(x,hanshu);
subplot(2,2,4);%绘制输出y的波形plot(y);
grid on;
axis([0,100,-2,2]);
xlabel('x');
ylabel('y');
title('输出信号y的波形')
function y=juanji(x,h)
lx=length(x);
lh=length(h);
sum=0;
lz=lx+lh;
for n=1:lz
for k=1:(lz-1)
if k<n&&(n-k)<=lh
if k<=lx
ji=x(k)*h(n-k);
sum=sum+ji;
end
end
end
if n>=2
y(n-1)=sum; sum=0;
end
end。