필터 지우기
필터 지우기

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

조회 수: 3 (최근 30일)
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
Sam Chak 2024년 6월 11일
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')

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider
Star Strider 2024년 6월 11일
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>
<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.
.

추가 답변 (1개)

sai charan sampara
sai charan sampara 2024년 6월 11일
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
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Vibration Analysis에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by