用双线性变换法设计一个数字巴特沃斯型带阻滤波器 用MATLAB编程 急！！！

%当1KHz≤f≤2KHz时，衰减至少为18 dB，当f≤500Hz，f≥3KHz时，波纹在3dB内，采样频率fs为10KHz。

%要求采用巴特沃斯型模拟滤波器进行设计，画出所设计滤波器的幅度特性曲线，并写出其系统函数H(z)的表达式。 f和pi之间怎么转化呢？f≤500Hz，f≥3KHz这在程序中怎么表示？求完整程序 急！！！

clear all

close all

clc

fs=10*1e+3

fsl=1*1e+3

fsh=2*1e+3

% fpl=500

fpl=800

fph=3*1e+3

alphal=18

alphah=3

ws1=fsl/fs

ws2=fsh/fs

wp1=fpl/fs

wp2=fph/fs

wp=[wp1 wp2]

ws=[ws1 ws2]

[n,Wn] = buttord(wp,ws,alphah,alphal)

[z,p,k] = buttap(n)

b=poly(z)

a=poly(p)

Wo=sqrt(2*pi*fpl*2*pi*fph)

Bw=2*pi*fph-2*pi*fpl

[bt,at] = lp2bs(b,a,Wo,Bw)

[numd,dend] = bilinear(bt,at,fs)

[H,w]=freqz(numd,dend)

plot(w/pi*fs/2,abs(H))

grid

说明：（1）为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽量使滤波器过渡带宽尽可能宽。

（2）与信号产生函数mstg相同，采样频率Fs=10kHz。

（3）为了滤波器阶数最低，选用椭圆滤波器。（之后，再依次实现巴特沃斯、切比雪夫1、切比雪夫2数字滤波器）

按照图2 所示的程序框图编写的实验程序为exp1.m。

2、实验程序清单

%实验1程序exp1.m

% IIR数字滤波器设计及软件实现

clear allclose all

Fs=10000T=1/Fs %采样频率

%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st

st=mstg

%低通滤波器设计与实现=========================================

fp=280fs=450

wp=2*fp/Fsws=2*fs/Fsrp=0.1rs=60 %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs)%调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp) %调用ellip计算椭圆带通DF系统函数系数向量B和A

y1t=filter(B,A,st)%滤波器软件实现

% 低通滤波器设计与实现绘图部分

figure(5)

subplot(2,1,1)

myplot(B,A) %调用绘图函数myplot绘制损耗函数曲线

yt='y_1(t)'

subplot(2,1,2)

tplot(y1t,T,yt)%调用绘图函数tplot绘制滤波器输出波形

%带通滤波器设计与实现====================================================

fpl=440fpu=560fsl=275fsu=900

wp=[2*fpl/Fs,2*fpu/Fs]ws=[2*fsl/Fs,2*fsu/Fs]rp=0.1rs=60

[N,wp]=ellipord(wp,ws,rp,rs) %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp)%调用ellip计算椭圆带通DF系统函数系数向量B和A

y2t=filter(B,A,st)%滤波器软件实现

% 带通滤波器设计与实现绘图部分

figure(3)

subplot(2,1,1)

myplot(B,A) %调用绘图函数myplot绘制损耗函数曲线

yt='y_2(t)'

subplot(2,1,2)

tplot(y2t,T,yt)%调用绘图函数tplot绘制滤波器输出波形

%高通滤波器设计与实现================================================

fp=890fs=600

wp=2*fp/Fsws=2*fs/Fsrp=0.1rs=60 %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs) %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp,'high')%调用ellip计算椭圆带通DF系统函数系数向量B和A

y3t=filter(B,A,st)%滤波器软件实现

% 高低通滤波器设计与实现绘图部分

figure(4)

subplot(2,1,1)

myplot(B,A) %调用绘图函数myplot绘制损耗函数曲线

yt='y_3(t)'

subplot(2,1,2)

tplot(y3t,T,yt)%调用绘图函数tplot绘制滤波器输出波形

function myplot(B,A)

%时域离散系统损耗函数绘图

%B为系统函数分子多项式系数向量

%A为系统函数分母多项式系数向量

[H,W]=freqz(B,A,1000)

m=abs(H)

plot(W/pi,20*log10(m/max(m)))grid on

xlabel('\omega/\pi')ylabel('幅度(dB)')

axis([0,1,-80,5])title('损耗函数曲线')

function tplot(xn,T,yn)

%时域序列连续曲线绘图函数

% xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串）

% T为采样间隔

n=0:length(xn)-1t=n*T

plot(t,xn)

xlabel('t/s')ylabel(yn)

axis([0,t(end),min(xn),1.2*max(xn)])

function st=mstg

N=2000

Fs=10000T=1/FsTp=N*T

t=0:T:(N-1)*Tk=0:N-1f=k/Tp

fc1=Fs/10

fm1=fc1/10

fc2=Fs/20

fm2=fc2/10

fc3=Fs/40

fm3=fc3/10

xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t)

xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t)

xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t)

st=xt1+xt2+xt3

fxt=fft(st,N)

subplot(3,1,1)

plot(t,st)gridxlabel('t/s')ylabel('s(t)')

axis([0,Tp/8,min(st),max(st)])title('(a) s(t)的波形')

subplot(3,1,2)

stem(f,abs(fxt)/max(abs(fxt)),'.')gridtitle('(b) s(t)的频谱')

axis([0,Fs/5,0,1.2])

xlabel('f/Hz')ylabel('幅度')

---