How to correct the recording of a daamaged accelerometer in earthquake analysis
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
charbel
il 11 Giu 2024
Commentato: Star Strider
il 11 Giu 2024
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 Commenti
Sam Chak
il 11 Giu 2024
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')
Risposta accettata
Star Strider
il 11 Giu 2024
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')
acceleration
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
Fstd = std(diff(acceleration(:,1))) % Sampling Frequency Variation
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 Commenti
Più risposte (1)
sai charan sampara
il 11 Giu 2024
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 Commenti
Image Analyst
il 11 Giu 2024
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.
Vedere anche
Categorie
Scopri di più su Earthquake Engineering in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!