How to correct the recording of a daamaged accelerometer in earthquake analysis

Question

charbel 2024-6-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2127361-how-to-correct-the-recording-of-a-daamaged-accelerometer-in-earthquake-analysis

评论： Star Strider 2024-6-11

acceleration.mat

Hello everyone,

I have the recording of two accelerometers but it turned out that one of them did not work properly and the recording is messed up (See photo, The two recordings should be the same). However I think if processed well I can recover the lost acclerogram. Does anyone have an idea how I can proceed. I have attached the recordings.

Thank you in advance

2 个评论
显示无隐藏无

Sam Chak 2024-6-11

在 MATLAB Online 中打开

acceleration.mat

Hi @charbel

This is an attempt to plot the results. As I do not have expertise in seismology, I cannot definitively determine which output is preferable and which is less so. However, a typical accelerogram would be expected to resemble the first (blue) trace.

load('acceleration.mat');

t = acceleration(:,1); % time

dat1= acceleration(:,2); % blue data

dat2= acceleration(:,3); % red data

subplot(211)

plot(t, dat1, 'color', "#0072BD"), grid on

xlim([t(1), t(end)]), ylim([-1, 2]), xlabel('Time'), ylabel('Acc / [cm/s^{2}]')

title('Accelerogram 1')

subplot(212)

plot(t, dat2, 'color', "#D95319"), grid on

xlim([t(1), t(end)]), ylim([-1, 2]), xlabel('Time'), ylabel('Acc / [cm/s^{2}]')

title('Accelerogram 2')

charbel 2024-6-11

Thank you for your time to answer

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2024-6-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2127361-how-to-correct-the-recording-of-a-daamaged-accelerometer-in-earthquake-analysis#answer_1470391

在 MATLAB Online 中打开

acceleration.mat

If you have the Signal Ptocessing Toolbox, an alternative approach would be to use the highpass (or bandpass to also eliminiate the high-frequency noise spike at about 16.5 Hz) function to filter out the low-frequency components from the ‘broken’ accerometer.

That would be something like this —

load('acceleration.mat')

whos('file','acceleration')

Name Size Bytes Class Attributes acceleration 24000x3 576000 double

acceleration

acceleration = 24000x3

0 0.0861 0.0806 0.0052 0.0868 0.0834 0.0104 0.0687 0.0651 0.0156 0.0393 0.0353 0.0208 0.0038 -0.0016 0.0260 -0.0243 -0.0321 0.0312 -0.0427 -0.0498 0.0365 -0.0429 -0.0493 0.0417 -0.0240 -0.0272 0.0469 0.0044 0.0015

<mw-icon class=""></mw-icon>

figure

plot(acceleration(:,1), acceleration(:,[2 3]))

grid

xlabel('Time')

ylabel('Amplitude')

title('Original Time Domain')

legend('Acceleromentr_1','Accelerometer_2', 'Location','best')

[FTs1,Fv] = FFT1(acceleration(:,[2 3]), acceleration(:,1));

figure

plot(Fv, abs(FTs1)*2)

grid

% set(gca, 'XScale','log')

xlabel('Frequency')

ylabel('Magnitude')

title('Fourier Transform')

xlim([0 20])

figure

plot(Fv, abs(FTs1)*2)

grid

% set(gca, 'XScale','log')

xlabel('Frequency')

ylabel('Magnitude')

title('Fourier Transform (Zoomed)')

xlim([0 1])

xline(0.45, '--k', 'Highpass Filter: F_{co} = 0.5 Hz')

Fs = 1/mean(diff(acceleration(:,1))) % Sampling Frequency

Fs = 192

Fstd = std(diff(acceleration(:,1))) % Sampling Frequency Variation

Fstd = 5.0181e-15

acc_hp_filt = highpass(acceleration(:,[2 3]), 0.5, Fs, 'ImpulseResponse','iir');

figure

plot(acceleration(:,1), acc_hp_filt)

grid

xlabel('Time')

ylabel('Amplitude')

title('Highpass Filtered Time Domain')

figure

plot(acc_hp_filt(:,1), acc_hp_filt(:,2))

grid

xlabel('Accelerometer_1 (Highpass Filtered)')

ylabel('Accelerometer_2 (Highpass Filtered)')

axis('equal')

acc_bp_filt = bandpass(acceleration(:,[2 3]), [0.5 7.5], Fs, 'ImpulseResponse','iir');

figure

plot(acceleration(:,1), acc_bp_filt)

grid

xlabel('Time')

ylabel('Amplitude')

title('Bandpass Filtered Time Domain')

figure

plot(acc_bp_filt(:,1), acc_bp_filt(:,2))

grid

xlabel('Accelerometer_1 (Bandpass Filtered)')

ylabel('Accelerometer_2 (Bandpass Filtered)')

axis('equal')

function [FTs1,Fv] = FFT1(s,t)

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

FTs1 = FTs(Iv,:);

end

It is necessary to use the Fourier transform in order to design the filtter passband frequencies correctly. This is overly detailed to show the relative effects of the two filter approaches. I filtered both signals with the same filters.

.

2 个评论
显示无隐藏无

charbel 2024-6-11

Thank you very much for this great response

Star Strider 2024-6-11

As always, my pleasure!

请先登录，再进行评论。

Answer 2

sai charan sampara 2024-6-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2127361-how-to-correct-the-recording-of-a-daamaged-accelerometer-in-earthquake-analysis#answer_1470181

在 MATLAB Online 中打开

acceleration.mat

Hello,

A possible solution can be to remove the mean over a small interval from the existing data to change the data to be centered around zero or around the mean of the correct data. It can be done similar to the code shown below:

load("acceleration.mat");

plot(acceleration(:,1),acceleration(:,2));

hold on

plot(acceleration(:,1),acceleration(:,3));

hold off

n=acceleration(:,3);

for i=1:10:length(acceleration(:,3))-9

n(i:i+9)=acceleration(i:i+9,3)-mean(acceleration(i:i+9,3))+mean(acceleration(:,2));

end

plot(acceleration(:,1),acceleration(:,2));

hold on

plot(acceleration(:,1),n,'g');

hold off

legend(["Accelogram 1","Corrected Accelogram 2"]);

The number of datapoints to be taken at a time for "mean" can be changed based on the requirement.

2 个评论
显示无隐藏无

Image Analyst 2024-6-11

To avoid changing the (probably) "good" data before the sensor went bad, you can start the correction just where it starts to go bad, like around x=50 or so, instead of at the first index.

charbel 2024-6-11