Fixing post load for acceleration signal with Savitzky-Golay filter

조회 수: 18 (최근 30일)
Emily Keys
Emily Keys 2023년 4월 7일
편집: E. Cheynet 2023년 4월 14일
Hi I have attached a code that processes acceleration data and applies a savitzky-Golay filter from another post however the post-load data is coming out incorrect as it slopes downwards where it should follow the orange line. The orange line is the accurate displacement measured from an LVDT and I am trying to match it using the trice integrated acceleration as can be seen in the photo below. My code is also attached. Thank you.
clear;clc;
Sheets = ["Sheet1","Sheet2","Sheet3","Sheet4","Sheet5","Sheet6","Sheet7","Sheet8","Sheet9"];
for test = 1 : 9
%% Set up the Import Options and import the data
opts = spreadsheetImportOptions("NumVariables", 1);
% Specify sheet and range
opts.Sheet = Sheets(test);
opts.DataRange = "C2:C246";
% Specify column names and types
opts.VariableNames = "VarName3";
opts.VariableTypes = "double";
% Import the data
tbl = readtable("C:\Users\Emily\LVDTData.xls", opts, "UseExcel", false);
%% Convert to output type
LDVT(:,test) = (59-(tbl.VarName3));
%% Clear temporary variables
clear opts tbl
end
clear test Sheets
fsL = 8.2;
tL = 0:(1/fsL):(length(LDVT(:,1))-1)/fsL;
tL = tL+5.06;
%% Now to unpack the JA and Kistler
%% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 4, "Encoding", "UTF-8");
% Specify range and delimiter
opts.DataLines = [5, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["Timestamp", "VarName2", "Timestamp1", "VarName4"];
opts.VariableTypes = ["double", "double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
tbl = readtable("C:\Users\Emily\Accelerations.csv", opts);
%% Convert to output type
tAcc = tbl.Timestamp;
KistlerWhole = tbl.VarName2;
JAWhole = tbl.VarName4;
fsA = 2048;
%% Clear temporary variables
clear opts tbl
AccStart = [ 65 ; 315; 655; 780; 990; 1399; 1502; 1600; 1844];
AccEnd = [ 105 ; 355; 695; 820; 1030; 1439; 1542; 1640; 1884];
for test = 1 : 9
[~,a] = min(abs(tAcc-AccStart(test)));
[~,b] = min(abs(tAcc-AccEnd(test)));
JA(:,test) = JAWhole(a:b);
Kistler(:,test) = KistlerWhole(a:b);
% JA(:,test) = detrend(JA(:,test));
% Kistler(:,test) = detrend(Kistler(:,test));
tA(:,test) = tAcc(a:b);
tA(:,test) = tA(:,test) -( tA(1,test) );
end
clear KistlerWhole JAWhole tAcc AccStart AccEnd test a b
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
hFig1 = figure;
tA = tA(:,9);
JA = JA(:,9);
plot(tA, JA, 'b.-', 'MarkerSize', 1);
grid on;
hold on;
fontSize = 20;
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
title('Original Signal', 'FontSize', fontSize);
hFig1.WindowState = 'maximized'; % Maximize the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);
% A moving trend is influenced by the huge outliers, so get rid of those first.
% Find outliers
outlierIndexes = isoutlier(JA);
plot(tA(outlierIndexes), JA(outlierIndexes), 'ro', 'MarkerSize', 15);
% Extract the good data.
tGood = tA(~outlierIndexes);
accelGood = JA(~outlierIndexes);
% plot(t(~outlierIndexes), accel(~outlierIndexes), 'mo', 'MarkerSize', 10); % Plot circles around the good data.
% Do a Savitzky-Golay filter (moving quadratic).
windowWidth = 51; % Smaller for tighter following of original data, bigger for smoother curve.
smoothedy = sgolayfilt(accelGood, 2, windowWidth);
hold on;
plot(tGood, smoothedy, 'r-', 'LineWidth', 2);
%legend('Original Signal', 'X axis', 'Outliers', 'Smoothed Signal');
% fill in the missing points.
smoothedy = interp1(tGood, smoothedy, tA);
% Now subtract the smoothed signal to get the variation
signal = JA - smoothedy;
% Plot it.
hFig2 = figure;
plot(tA, signal, 'b.-', 'MarkerSize', 9);
grid on;
hold on;
title('Corrected Signal', 'FontSize', fontSize);
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
hFig2.Units = 'normalized';
hFig2.Position = [.2, .2, .5, .5]; % Size the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);
VelCor = cumtrapz(tA,signal);
DispCorr = cumtrapz(tA,VelCor);
DispCorr = DispCorr*10^4;
figure(); clf; hold on
subplot(1,1,1);
plot(tA, DispCorr, 'b');hold on
plot(tL,LDVT(:,9));hold off
grid on;
hold on;
title('Corrected Displacement', 'FontSize', fontSize);
xlabel('Time', 'FontSize', fontSize);
ylabel('Displacement', 'FontSize', fontSize);
hFig2.Units = 'normalized';
hFig2.Position = [.2, .2, .5, .5]; % Size the figure window.
% Draw a line at y=0
  댓글 수: 3
Emily Keys
Emily Keys 2023년 4월 7일
Hi Mathieu, Yes that code is very helpful so thank you for that, but for my work I am trying out multiple ways to fix the issues of drift and low frequency noise as to compare.
Mathieu NOE
Mathieu NOE 2023년 4월 11일
ok
I'll be interested to see what kind of alternatives can work, because so far , all posts I have seen or answered on this topic are using the same approach (more or less)

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

답변 (1개)

E. Cheynet
E. Cheynet 2023년 4월 14일
편집: E. Cheynet 2023년 4월 14일
There are essentially two main approaches to converting acceleration to displacement signal: (1) double integration as wisely mentioned by Mathieu Noe or (2) using the fast Fourier transform. Here is an example in Matlab File Exchange. I am not sure which method works best. At a point, it is simply no longer possible to retrieve the quasi-static displacement from the accelerometer records, unless a better (and more expensive) sensor is used.

카테고리

Help CenterFile Exchange에서 Matched Filter and Ambiguity Function에 대해 자세히 알아보기

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by