%Samp10_3
load hns.dat ;                           %假設的地震儀記錄之前的地面運動
Xt=hns;                                 %把數據賦值給變量
Fs=50;                                   %設定采樣率 單位(HZ)
dt=1/Fs;                                 %求采樣間隔 單位(s)
N=length(Xt);                            %得到序列的長度
t=[0:N-1]*dt;                            %時間序列
%設計一個橢圓濾波器，假定為寬頻帶地震儀
ws=[0.00001 24.9]*2/Fs;wp=[0.001 24.8]*2/Fs;   %通帶和阻帶邊界頻率（歸一化頻率）
Rp=1;Rs=50;Nn=512; %通帶波紋和阻帶衰減以及繪制頻率特性的數據點數
[Order,Wn]=ellipord(wp,ws,Rp,Rs);  %求取數字濾波器的最小階數和歸一化截止頻率
[b,a]=ellip(Order,Rp,Rs,Wn); %按最小階數、截止頻率、通帶波紋和阻帶衰減設計濾波器
figure(1)
[H,f]=freqz(b,a,Nn,Fs); %按傳遞函數系數、數據點數和采樣頻率求得濾波器的頻率特性
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;
y=filtfilt(b,a,Xt);  %在寬帶濾波器上的輸出
figure(2)
subplot(2,1,1),plot(t,Xt)
xlabel('時間/s'),title('輸入信號');
ylabel('振幅');
subplot(2,1,2),plot(t,y)
xlabel('時間/s'),title('輸出信號');
ylabel('振幅');
%已知寬頻帶地震儀的頻率特性,恢復地面運動
[H,f]=freqz(b,a,N,Fs,'whole');  %得到地震儀的特性
Y=fft(y);
XX=zeros(1,N);
for ii=1:N      %按（10-7）式得到地面運動的頻率域表示
    if (H(ii)>1.0e-4)
        XX(ii)=Y(ii)./H(ii);
    end
end
disp=real(ifft(XX));   %得到地面運動
figure(3),plot(t,Xt,t,disp,'r')
legend('原始信號','恢復地面運動',1)
xlabel('時間/s')
ylabel('振幅');
%仿真到短周期地震儀上，短周期地震儀用一個窄帶橢圓濾波器來表示
ws=[0.01 4.5]*2/Fs;wp=[0.1 3.8]*2/Fs;   %通帶和阻帶邊界頻率（歸一化頻率）
Rp=1;Rs=20;Nn=512; %通帶波紋和阻帶衰減以及繪制頻率特性的數據點數
[order,Wn]=ellipord(wp,ws,Rp,Rs);  %求取數字濾波器的最小階數和歸一化截止頻率
[b,a]=ellip(order,Rp,Rs,Wn); %按最小階數、截止頻率、通帶波紋和阻帶衰減設計濾波器
figure(4)
[H1,f]=freqz(b,a,Nn,Fs); %按傳遞函數系數、數據點數和采樣頻率求得濾波器的頻率特性
subplot(2,1,1),plot(f,20*log10(abs(H1)))
xlabel('頻率/Hz');ylabel('振幅/dB');grid on;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H1)))
xlabel('頻率/Hz');ylabel('相位/^o');grid on;
figure(5)
y1=filtfilt(b,a,Xt);  %在窄帶濾波器上的輸出
[H1,f]=freqz(b,a,N,Fs,'whole');  %得到地震儀的特性
XX1=zeros(1,N);
for ii=1:N    %按（10-6）式得到仿真結果
if (abs(H(ii))>1.0e-4)
    XX1(ii)=Y(ii).*H1(ii)/H(ii);
end
end
x1=ifft(XX1);
plot(t,y1,t,real(x1),'r') %繪制輸入信號
legend('實際輸出','仿真輸出',1)
xlabel('時間/s');
ylabel('振幅');
grid on;
%仿真到長周期地震儀上，長周期地震儀用一個窄帶橢圓濾波器來表示
ws=0.1*2/Fs;wp=0.02*2/Fs;   %通帶和阻帶邊界頻率（歸一化頻率）
Rp=1;Rs=30;Nn=512; %通帶波紋和阻帶衰減以及繪制頻率特性的數據點數
[Order,Wn]=ellipord(wp,ws,Rp,Rs);  %求取數字濾波器的最小階數和歸一化截止頻率
[b,a]=ellip(Order,Rp,Rs,Wn); %按最小階數、截止頻率、通帶波紋和阻帶衰減設計濾波器
figure(6)
y1=filtfilt(b,a,Xt);  %在窄帶濾波器上的輸出
[H1,f]=freqz(b,a,N,Fs,'whole');  %得到地震儀的特性
XX1=zeros(1,N);
for ii=1:N
if (abs(H(ii))>1.0e-4)   %為了防止H值太小將該頻率的信號放大
    XX1(ii)=Y(ii).*H1(ii)./H(ii);  %按（10-6）式得到仿真結果
end
end
x1=ifft(XX1);
plot(t,y1,t,real(x1),'r:') %繪制輸入信號
legend('實際輸出','仿真輸出',1)
xlabel('時間/s');
ylabel('振幅');
grid on;