%Appl6_1
load ml031025ns.txt   %加載地震波形記錄
dt=0.02;   %中國數字地震臺網的采樣間隔為0.02s,采樣頻率為50Hz
x=ml031025ns';
wp=0.5*2*dt;  ws=0.2*2*dt; %根據采樣頻率將濾波器邊界頻率進行轉換
Rp=1;Rs=30;    %通帶波紋和阻帶衰減
Nn=128;      %顯示濾波器頻率特性的數據長度
[N,Wn]=buttord(wp,ws,Rp,Rs);  %求得濾波器的最小階數和截止頻率（歸一化頻率）
[b,a]=butter(N,Wn,'high');  %設計Butterworth高通數字濾波器
figure(1)
[H,f]=freqz(b,a,Nn,1/dt);    %用Nn點求得頻率特性
subplot(2,1,1),plot(f,20*log10(abs(H)));
xlabel('頻率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('頻率/Hz');ylabel('相位/^o');grid on;
n=0:length(x)-1;
t=n*dt;   %時間序列
figure(2)
subplot(2,1,1),plot(t,x); title('輸入信號')  %繪制輸入信號
y=filter(b,a,x);   %對輸入信號進行濾波
subplot(2,1,2),plot(t,y),title('輸出信號') %繪制輸出信號
xlabel('時間/s')