%用ellip函數設計濾波器
[y,fs,bits]=wavread('d:\test.wav');  %將文件讀進來
%從這里開始設計濾波器
wp=1000/(fs/2);ws=1200/(fs/2);
[N,wc]=ellipord(wp,ws,1,100);
[b,a]=ellip(N,1,100,wc);
%[bz,az]=bilinear(b,a,fs); %映射為數字的
figure(1);
freqz(b,a,512,fs/2); 
figure(2);
subplot(2,3,1);
plot(y,'LineWidth',2);
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) ;
xlabel('k')                     %把x軸定義為w
title('原信號幅度')                   %標明此圖形畫的是幅度
subplot(2,3,3)
plot(t2,spimag)  ;  
xlabel('k')                     %把x軸定義為虛部
title('原信號相位')                   %把x軸定義為虛部
y0=filter(b,a,y);
n=length(y0);               %求得經過低通濾波器的y0的長度
xx=0:n-1;
subplot(2,3,4);
plot(xx,y0);
xlabel('n')
ylabel('幅度');
title('經過低通濾波器后信號幅度')
y1=fft(y0);
y0am=abs(y1);
subplot(2,3,5);
plot(xx,y0am);
title('經過低通濾波器后信號幅度譜');
subplot(2,3,6);
yimag=angle(y1);
plot(xx,yimag);
xlabel('n');
ylabel('相位');
title('經過低通濾波器后信號相位')
sound(y0,fs)