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])