%Appl7_6
load bssud.txt   %加載數據
x=bssud; %輸入信號
n=320;    %濾波器階數
dt=0.02;   %中國數字地震臺網的采樣間隔為0.02s,采樣頻率為50Hz
t=[0:length(x)-1]*dt;    %時間序列
f=[0 0.5/25 0.6/25 1.5/25 1.6/25 1];            %頻率向量
a=[0 0 1 1 0 0];                 %振幅向量
b=firls(n,f,a);              %采用 firls 設計濾波器
[h,w]=freqz(b);                %計算其頻率響應
bb=remez(n,f,a);             %采用 remez 設計濾波器
[hh,w]=freqz(bb);        
figure(1)
subplot (2,1,1),plot (t,x),title ('輸入信號')           %繪出輸入信號波形
subplot (2,1,2),plot(w/pi,abs(h),w/pi,abs(hh),'r:');         %繪制濾波器的幅頻響應
xlabel('歸一化頻率');ylabel('振幅');
legend('firls','remez');                       %給出圖例
%濾波
y1=filtfilt(b,1,x); %產生零相位輸出，采用filtfilt 對輸入信號濾波
y2=filtfilt(bb,1,x); %產生零相位輸出，采用filtfilt 對輸入信號濾波
t=[0:length(x)-1]*dt;
figure(2)
subplot (2,1,1),plot (t,y1,'b')                  %繪出輸出信號波形
title('輸出信號'),xlabel('時間/s')
legend('firls');                                 %給出圖例
subplot (2,1,2),plot (t,y2,'r')                  %繪出輸出信號波形
title('輸出信號'),xlabel('時間/s')
legend('remez');                                 %給出圖例