Fixing post load for acceleration signal with Savitzky-Golay filter
4 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
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
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.
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!