kuandaizhaidai.m

用真实的地震数据和数字滤波器方法来揭示地面运动向窄带上的仿真。

%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;