高分求：用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器

我当时改写的课程设计程序，希望对你有帮助（参数自己改，很容易的）
%%%%%%%%%%%%% 低通滤波 %%%%%%%%%%%%%%%%
clear;
clear clf;
%%% 对连续时间信号进行采样
f1=2;f2=5;f3=8;
fs=20;Ts=1/fs;
M=200;
k=0:M-1;
fk=cos(2*pi*f1*k*Ts)+cos(2*pi*f2*k*Ts)+cos(2*pi*f3*k*Ts);
%figure(1)
subplot(411)
plot(k,fk)%stem(k,fk)
xlabel ' '
title '滤波前的波形图';

N = M;
F = fft(fk, N);
subplot(412)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(F(1:N/2))/N);
xlabel ' '
title '滤波前的频谱曲线';

h=[0.00111829516864 -0.00389476479172 -0.01603491745519 -0.02036377118215 0.02095180705130 0.12449781344246...
0.24450683184615 0.29843741184102 0.24450683184615 0.12449781344246 0.02095180705130 -0.02036377118215...
-0.01603491745519 -0.00389476479172 0.00111829516864];
yk = conv(fk,h);

%figure(2)
subplot(413)
plot(0:M+15-2,yk,'g')%stem(0:M+15-2,yk)
xlabel ' '
title '低通滤波后的波形图';
axis([0 M -1 1])
Y = fft(yk, N);
subplot(414)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(Y(1:N/2))/N,'g');
title '低通滤波后的频谱曲线';
%===================注：与高通滤波不同之处在于h的取值
%%%%%%%%%%%%% 课程设计(2) %%%%%%%%%%%%%
%%%%%%%%%%%%% 高通滤波 %%%%%%%%%%%%%%%%
clear;
clear clf;
%%% 对连续时间信号进行采样
f1=2;f2=5;f3=8;
fs=20;Ts=1/fs;
M=200;
k=0:M-1;
fk=cos(2*pi*f1*k*Ts)+cos(2*pi*f2*k*Ts)+cos(2*pi*f3*k*Ts);
%figure(1)
subplot(411)
plot(k,fk)
xlabel ' '
title '滤波前的波形图';

N = M;
F = fft(fk, N);
subplot(412)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(F(1:N/2))/N);
xlabel ' '
title '滤波前的频谱曲线';

h=[-0.00111829516864 -0.00389476479172 0.01603491745519 -0.02036377118215 -0.02095180705130 0.12449781344246...
-0.24450683184615 0.29843741184102 -0.24450683184615 0.12449781344246 -0.02095180705130 -0.02036377118215...
0.01603491745519 -0.00389476479172 -0.00111829516864];
yk = conv(fk,h);

%figure(2)
subplot(413)
plot(0:M+15-2,yk)
xlabel ' '
title '滤波后的波形图';
axis([0 M -1 1])
Y = fft(yk, N);
subplot(414)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(Y(1:N/2))/N);
title '滤波前的频谱曲线';
%%%%%%%%%%%%% 课程设计(4) %%%%%%%%%%%%%
%%%%%%%%%%%%% 带通滤波 %%%%%%%%%%%%%%%%
clear;
clear clf;
%%% 对连续时间信号进行采样
f1=2;f2=5;f3=8;
fs=20;Ts=1/fs;
M=200;
k=0:M-1;
fk=cos(2*pi*f1*k*Ts)+cos(2*pi*f2*k*Ts)+cos(2*pi*f3*k*Ts);
subplot(411)
plot(k,fk)
xlabel ' '
title '滤波前的波形图'

N = M;
F = fft(fk, N);
subplot(412)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(F(1:N/2))/N);
xlabel ' '
title '滤波前的频谱曲线'

h=[ 0 0.00809904403983 0 0.04234583818052 0 -0.25888938815435 0 0.41372763540994 0 -0.25888938815435 0 0.04234583818052 0 0.00809904403983 0];
yk = conv(fk,h);

figure(2)
subplot(413)
plot(0:M+15-2,yk)
axis([0 M -1 1])
Y = fft(yk, N);
subplot(414)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(Y(1:N/2))/N);
title '滤波前的频谱曲线'
%%%%%%%%%%%%% 课程设计(4) %%%%%%%%%%%%%
%%%%%%%%%%%%% 带阻滤波 %%%%%%%%%%%%%%%%
clear;
clear clf;
%%% 对连续时间信号进行采样
f1=2;f2=5;f3=8;
fs=20;Ts=1/fs;
M=200;
k=0:M-1;
fk=cos(2*pi*f1*k*Ts)+cos(2*pi*f2*k*Ts)+cos(2*pi*f3*k*Ts);

subplot(411)
plot(k,fk)
xlabel ' '
title '滤波前的波形图'
N = M;
F = fft(fk, N);
subplot(412)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(F(1:N/2))/N);
xlabel ' '
title '滤波前的频谱曲线'

h=[ 0 -0.00780645449547 0 -0.04081603423850 0 0.24953663889817 0 0.59817169967160...
0 0.24953663889817 0 -0.04081603423850 0 -0.00780645449547 0];
yk = conv(fk,h);

subplot(413)
plot(0:M+15-2,yk,'r')
xlabel ' '
title '带阻滤波后前的波形图'
axis([0 M -1 1])
Y = fft(yk, N);
subplot(414)
plot(2*pi*(0:N/2-1)/N/pi, 2*abs(Y(1:N/2))/N,'r');
title '带阻滤波前的频谱曲线'

双线性变换法设计IIR带通滤波器，其中带通的中心频率为ωp0=0.5π，通带截止频率ωp1=0.4π,ωp2=0.6π;通带最大衰减αp=3dB;阻带最小衰减αs=15dB;阻带截止频率ωs2=0.7π,程序如下：

clear
wp0=0.5*pi;wp1=0.4*pi;wp2=0.6*pi;
Ap=3;ws2=0.7*pi;As=15;T=2; %数字带通滤波器技术指标
ws1=wp0-(ws2-wp0); %计算带通滤波器的阻带下截止频率
wc1=(2/T)*tan(wp1/2);wc2=(2/T)*tan(wp2/2);
wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);
w0=(2/T)*tan(wp0/2); %频率预畸变
B=wc2-wc1; %带通滤波器的通带宽度
normwr1=(((wr1^2)-(w0^2))/(B*wr1));
normwr2=(((wr2^2)-(w0^2))/(B*wr2));
normwc1=(((wc1^2)-(w0^2))/(B*wc1));
normwc2=(((wc2^2)-(w0^2))/(B*wc2)); %带通到低通的频率变换
if abs(normwr1)>abs(normwr2)
normwr=abs(normwr2)
else normwr=abs(normwr1)
end
normwc=1; %将指标转换成归一化模拟低通滤波器的指标
[N,wn]=ellipord(normwc,normwr,Ap,As,'s');%设计归一化的模拟低通滤波器阶数N和3db截止频率
[bLP,aLP]=butter(N,normwc,normwr,'s'); %计算相应的模拟滤波器系统函数G(p)
[bBP,aBP]=lp2bp(bLP,aLP,w0,B); %模拟域频率变换，将G(P)变换成模拟带通滤波器H(s)
[b,a]=bilinear(bBP,aBP,1/T); %用双线性变换法将H(s)转换成数字带通滤波器H(z)
w=[0:500]*2*pi/500;
freqz(b,a,w);

wp=0.2*pi;ws=0.3*pi;
Fs=4000;T=1/Fs;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
rp=1;rs=15;as=15;
ripple=10^(-rp/20);attn=10^(-rs/20);
[n,wn]=buttord(OmegaP,OmegaS,rp,rs,'s');
[z,p,k]=Buttap(n);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2lp(b,a,wn);
[b,a]=bilinear(bt,at,Fs);
[db,mag,pha,grd,w]=freqz_m(b,a);
%
%下面绘出各条曲线
subplot(2,2,1);plot(w/pi,mag);title('Magnitude Frequency幅频特性');
xlabel('w(/pi)');ylabel('|H(jw)|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[0 attn ripple 1]);grid

subplot(2,2,2);plot(w/pi,db);title('Magnitude Frequency幅频特性(db)');
xlabel('w(/pi)');ylabel('dB');
axis([0,1,-30,5]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[-60 -as -rp 0]);grid

subplot(2,2,3);plot(w/pi,pha/pi);title('Phase Frequency相频特性');
xlabel('w(/pi)');ylabel('pha(/pi)');
axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);grid

subplot(2,2,4);plot(w/pi,grd);title('Group Delay群延时');
xlabel('w(/pi)');ylabel('Sample');
axis([0,1,0,15]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);grid

p1=400;
p2=500;%通带边界频率
s1=350;
s2=550;%阻带截止频率
Ap=1;%通带最大衰减
As=40;%阻带最小衰减
Ft=2000;%抽样频率
T=2;
wp1=2*pi*p1/Ft;
wp2=2*pi*p2/Ft;
ws1=2*pi*s1/Ft;
ws2=2*pi*s2/Ft;
Wp1=(2/T)*tan(wp1/2);
Wp2=(2/T)*tan(wp2/2);
Ws1=(2/T)*tan(ws1/2);
Ws2=(2/T)*tan(ws1/2)
W0=Wp1*Wp2;
w0=sqrt(W0);
BW=Wp2-Wp1; %带通滤波器的通带宽度
lp=1; %归一化处理
ls=Ws1*BW/(W0-Ws1^2);
[N,Wn]=ellipord(lp,ls,Ap,As,'s');
[B,A]=ellip(N,1,40,Wn,'s');
[BT,AT]=lp2bp(B,A,w0,BW);
[num,den]=bilinear(BT,AT,0.5);
[z,p,k]=tf2zp(num,den);
figure(1);
zplane(z,p);
title('零极点')
[h,w]=freqz(num,den,512);
figure(2)
plot(w/pi,20*log10(abs(h)));
axis([0 1 -100 1]);
title('频谱特性曲线')
grid
n=0:800;k=n/8000;%通过滤波器3
f1=2*pi*450;
f2=2*pi*6000;
x=sin(f1*k)+sin(f2*k);
y=filter(num,den,x);
x1=sin(f1*k)
figure(3)
plot(x1);%x1图形输出
axis([0,100*pi,-5,5]);
title('x1(t)');
x2=sin(f2*k);
figure(4)
plot(x2);%x2图形输出
axis([0,100*pi,-5,5]);
title('x2');
figure(5)
plot(x);
axis([0,100*pi,-5,5]);
title('输入信号');
figure(6)
plot(y);
axis([0,100*pi,-5,5]);
title('输出信号');

我试过了，可以的,数据改一下就行