%設(shè)計(jì)Butterworth低通濾波器
%數(shù)字濾波器指標(biāo)w=Omega.*T=Omega/fs：
[x1]=wavread('音樂(lè).wav')
wp=2*pi*1000/22050;     %數(shù)字通帶頻率(弧度)
ws=2*pi*1200/22050;     %數(shù)字阻帶頻率(弧度)
Rp=1;
Rs=100;
%預(yù)畸變：
fs=22050;T=1/fs;

OmegaP=2*tan(wp/2)/T;   %預(yù)修正原型通帶頻率
OmegaS=2*tan(ws/2)/T;   %預(yù)修正原型阻帶頻率
%計(jì)算模擬Butterworth濾波器：
[n,Wn]=buttord(OmegaP,OmegaS,Rp,15,'s');

[b,a]=butter(n,Wn,'s');
[bz,az]=bilinear(b,a,fs);%雙線性變換
[H,W]=freqz(bz,az,512);   %頻率響應(yīng)

%畫(huà)圖
figure(1)
title('FIR高通濾波器');
plot(W*fs/(2*pi),abs(H))
grid;
f2=filter(bz,az,x1);%將信號(hào)入濾波器

figure(2)
subplot(2,1,1)
plot(x1)
title('FIR低通濾波器濾波前的時(shí)域波形');
subplot(2,1,2)
plot(f2);
title('FIR低通濾波器濾波后的時(shí)域波形');
F0=fft(f2,1024);
f=fs*(0:511)/1024;

figure(3)
y2=fft(x1,1024);
subplot(2,1,1);
plot(f,abs(y2(1:512)));
title('FIR低通濾波器濾波前的頻譜')
xlabel('頻率/Hz');
ylabel('幅值');
subplot(2,1,2)
F2=plot(f,abs(F0(1:512)));
title('FIR低通濾波器濾波后的頻譜')
xlabel('頻率/Hz');
ylabel('幅值'); 
sound(f2,44100);    %濾波后的語(yǔ)音信號(hào)