%用fir1函數設計低通濾波器
 clf;
[y,fs,bits]=wavread('d:\test.wav');  %將文件讀進來
%設計低通濾波器
wp=1000/(fs/2);              %計算通帶的頻率
ws=1200/(fs/2);              %計算阻帶的頻率
fp=[1000 1200];              %通帶頻率
ap=10^0.05-1;                %求通帶衰減
as=10^(-5);                 %求阻帶衰減
dev=[ap as];
mags=[1 0];
[n,wn,beta,ftype]=kaiserord(fp,mags,dev,fs);         %求的用kaiserord時的濾波器階數
b=fir1(n,wn,ftype,kaiser(n+1,beta),'noscale');        %取濾波器的階數為n
figure(1);
freqz(b)   ;%畫出濾波器的幅度相應和相位
figure(2);
subplot(2,3,1);
plot(y,'LineWidth',2);
%axis([0 50000 -0.05 0.06]);
xlabel('時間/t')  ; ylabel('幅度/A');
title('原信號波形');   
spectrum=fft(y);                %對采樣點進行fft變換
n2=length(spectrum);             %計算fft變換后向量的長度
t2=0:n2-1;   
subplot(2,3,2);
spreal=abs(spectrum)   ;        %計算fft變換后序列的實部
spimag=angle(spectrum)  ;         %計算fft變換后序列的虛部
plot(t2,spreal)  ;
%axis([0 60000 -0.05 4]);  
xlabel('k')  ;                   %把x軸定義為w
title('原信號幅度');                   %標明此圖形畫的是幅度
subplot(2,3,3);
plot(t2,spimag)  ;  
xlabel('k')   ;                  %把x軸定義為虛部
title('原信號相位')     ;              %把x軸定義為虛部
y0=fftfilt(b,y);            %用低通濾波器對抽樣信號進行濾波
n=length(y0);
xx=0:n-1;
subplot(2,3,4);
plot(xx,y0);
%axis([0 50000 -0.05 0.06]);
xlabel('n');
ylabel('幅度');
y1=fft(y0);                   %對濾波器輸出信號進行fft變換
d=length(y1);
xx=0:d-1;
y0am=abs(y1);
subplot(2,3,5);
plot(xx,y0am);
%axis([0 60000 -0.05 4]); 
title('經過低通濾波器后信號幅度');
y0im=angle(y1);
subplot(2,3,6);
plot(xx,y0im);
axis([0 100000 -4 4]);
title('經過低通濾波器后信號相位');
sound(y0,fs);                            %回放錄制信號