x1=wavread('音樂.wav');
Rp=1;
Rs=8;
Fs=22050;
Ts=1/Fs;
wp1=2*1200*pi/Fs;ws1=2*1000*pi/Fs;
wp2=2*3000*pi/Fs;ws2=2*3200*pi/Fs;
wp=[wp1,wp2];
ws=[ws1,ws2];
wp11=2/Ts*tan(wp1/2); 
wp22=2/Ts*tan(wp2/2); 
BW=wp22-wp11; W0=wp11*wp22; W00=sqrt(W0); 
wpp=2/Ts*tan(wp/2);                 %將模擬指標(biāo)轉(zhuǎn)換成數(shù)字指標(biāo)
wss=2/Ts*tan(ws/2);
[N,Wn]=buttord(wpp,wss,Rp,Rs,'s');  %選擇濾波器的最小階數(shù)
[Z,P,K]=buttap(N);                  %創(chuàng)建巴特沃斯模擬濾波器
[Bap,Aap]=zp2tf(Z,P,K);        %變系統(tǒng)零點(diǎn)增益為傳遞函數(shù)形式
[b,a]=lp2bp(Bap,Aap,W00,BW);     %低通到帶通模擬濾波器的變換
[bz,az]=bilinear(b,a,Fs);           %用雙線性變換法實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換
[H,W]=freqz(bz,az);                 %繪制頻率響應(yīng)曲線


figure(1)
plot(W*Fs/(2*pi),abs(H))
grid
xlabel('頻率／Hz')
ylabel('頻率響應(yīng)幅度')
title('巴特沃斯帶通慮波器 ')

f1=filter(bz,az,x1);

figure(2)
subplot(2,1,1)
plot(x1)                          %畫出濾波前的時(shí)域圖
title('濾波前的時(shí)域波形');
subplot(2,1,2)
plot(f1);                         %畫出濾波后的時(shí)域圖
title('濾波后的時(shí)域波形');
sound(f1,44100);                    %播放濾波后的信號(hào)
F0=fft(f1,1024);
f=fs*(0:511)/1024;

figure(3)
y2=fft(x1,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512)));             %畫出濾波前的頻譜圖
title('濾波前的頻譜')
xlabel('Hz');
ylabel('幅值');
subplot(2,1,2)
F1=plot(f,abs(F0(1:512)));          %畫出濾波后的頻譜圖
title('濾波后的頻譜')
xlabel('Hz');
ylabel('幅值');