matlab 地震波过滤器设计

1、内容简介

略
550-可以交流、咨询、答疑

2、内容说明

略

3、仿真分析

close all
clear
clc
rng default
%% read data
dir_path = 'AOM0031811022001/';
filename = 'AOM0031811022001';
[ stream, stream_data ] = parseStreamFromData(dir_path,filename);
fs = stream.sampling_rate;
t = (0:length(stream_data{1,1})-1)/fs;
data1 = stream_data{1,1};
figure
plot(t, data1)
xlabel 时间/s
ylabel 幅值
title Origin
[f1,y1]=change_fft(fs,data1);
figure
plot(f1,y1)
xlabel 频率/Hz
ylabel 幅值
title Origin
%% define the noise
SNR = 15;
s = data1;

x=awgn(data1,SNR,'measured'); %观测信号 x=s+v.，给正弦波信号加入信噪比为-3dB的高斯白噪声
NOISE=x - s; % 高斯白噪声，误差信号，每次运行都不一样
figure
plot(t, NOISE)
xlabel 时间/s
ylabel 幅值
title noise

Y=data1+NOISE;
figure
plot(t, Y)
xlabel 时间/s
ylabel 幅值
title addNoise
[f1_noise,y1_noise]=change_fft(fs,Y);
figure
plot(f1_noise,y1_noise)
xlabel 频率/Hz
ylabel 幅值
title addNoise
%% use the kalman filter
Z=Y;
X=Z(1);
Q=0.1 ;
R=0.01;
P=0.01;
A=1;
H=1;
Xn=KaermanFilter(X,Z, A, Q, H, R,P);
figure
plot(t, data1, t, Y, t, Xn)
xlabel 时间/s
ylabel 幅值
legend('origin','before','after')
[f11,y11]=change_fft(fs,Xn);
figure
plot(f1,y1,f1_noise,y1_noise,f11,y11)
xlabel 频率/Hz
ylabel 幅值
legend('origin','before','after')
 

 

4、参考论文

略