High pass filter on cumtrapz

Hi! I am new to signal processing and in Matlab so my question may look dumb, but I was having trouble on the getting the graph that i want on my velocity and displacement. I want to use a highpass filter and not detrend code so i wont have to edit the order everytime i input another data. I attached the example file that i was working on too. I would really appreciate the help.

clear all;
close all;
clc;

%raw data from XYZ direction
NormChest = readtable('Chest-Normal.csv');                   %%%get the raw data 
C = table2array(NormChest(:,1));                               
x = table2array(NormChest(:,2));
y = table2array(NormChest(:,3));
z = table2array(NormChest(:,4));
time = C(:,1);
mag = table2array(NormChest(:,5));

%%design for low pass filter
fs = 1000; %sampling frequency                                 %%%trial first
fc = 20; %cut-off frequency                                    %%%trial first
order = 2; %6th order filter                                   %%%removes any low frequency

%%Filter acceleration signals
[b1, a1] = butter (order, fc/(fs/2));
acceleration  = filtfilt (b1, a1, mag);
figure (2)
plot (time, mag); hold on;
plot (time, acceleration, 'r'); hold on; axis tight;
xlabel('Time'); ylabel('Acceleration'); 
legend ('Raw', 'Filtered', 'Orientation','vertical', 'Location','northeast')

%first integration (acceleration-velocity)
fcv = 6;    %cut-off frequency
orderv = 6; %6th order filter
[b2, a2] = butter (orderv, fcv, "high");
fvtool = ([b2, a2]);
velocity = cumtrapz(time, acceleration);
plot (time, velocity);axis tight;
xlabel('Time'); ylabel('Velocity'); 

%%second integration (velocity-displacement)
displacement = cumtrapz(time, velocity);
%displacement = detrend (displacement, 6);                     %%%detrend slope graph
figure (4)
plot (time, displacement); axis tight;
xlabel('Time'); ylabel('Displacement');
time = time - time(1);

Answer 1

[Star Strider]

The ‘acceleration’ signal has a constant offset of about 1.03, so the integration is going to emphasize the offset and not the instantaneous acceleration. Subtracting that first:

mag = mag - mean(mag); % Remove D-C Offset

may get you closer to what you want.

Also, you do not use the highpass filter anywhere. That aside, I changed it so it will work:

[b2, a2] = butter (orderv, fcv*2/fs, "high"); % <— CHANGED

It probably does not matter here, however it is generally better to use the zero-pole-gain output ([z,p,k]) of butter and similar functions, and then use zp2sos to use a second-order-section implementation of the filter for stability. See that documentation for details.

Try this —

clear all;
close all;
clc;

%raw data from XYZ direction
NormChest = readtable('Chest-Normal.csv');                   %%%get the raw data 
C = table2array(NormChest(:,1));                               
x = table2array(NormChest(:,2));
y = table2array(NormChest(:,3));
z = table2array(NormChest(:,4));
time = C(:,1);
mag = table2array(NormChest(:,5));

mag = mag - mean(mag);                                          % Remove D-C Offset

%%design for low pass filter
fs = 1000; %sampling frequency                                 %%%trial first
fc = 20; %cut-off frequency                                    %%%trial first
order = 2; %6th order filter                                   %%%removes any low frequency

%%Filter acceleration signals
[b1, a1] = butter (order, fc/(fs/2));
acceleration  = filtfilt (b1, a1, mag);
figure (2)
plot (time, mag); hold on;
plot (time, acceleration, 'r'); hold on; axis tight;
xlabel('Time'); ylabel('Acceleration'); 
legend ('Raw', 'Filtered', 'Orientation','vertical', 'Location','northeast')

%first integration (acceleration-velocity)
fcv = 6;    %cut-off frequency
orderv = 6; %6th order filter
[b2, a2] = butter (orderv, fcv*2/fs, "high");                   % <— CHANGED
fvtool = ([b2, a2]);
velocity = cumtrapz(time, acceleration);

figure
plot (time, velocity);axis tight;
xlabel('Time'); ylabel('Velocity'); 

%%second integration (velocity-displacement)
displacement = cumtrapz(time, velocity);
%displacement = detrend (displacement, 6);                     %%%detrend slope graph
figure (4)
plot (time, displacement); axis tight;
xlabel('Time'); ylabel('Displacement');

time = time - time(1);

.

[Stella]

Thank you!! i was able to see the oscillation with these and used filter for a better graph, then also used the same to my displacement with different order. I really appreciate your fast response!!!

[Star Strider]

As always, my pleasure!

Answer 2

[Hassaan]

clear all;
close all;
clc;

% Load data from the CSV file
NormChest = readtable('Chest-Normal.csv');
C = table2array(NormChest(:,1));                               
x = table2array(NormChest(:,2));
y = table2array(NormChest(:,3));
z = table2array(NormChest(:,4));
time = C(:,1);
mag = table2array(NormChest(:,5));

% Design a high-pass filter to remove low-frequency drift
fs = 1000; % Sampling frequency
fchp = 0.5; % Cut-off frequency for high-pass filter, choose based on your data
order = 2; % Define the order of the filter, e.g., 1, 2, etc.
[b1, a1] = butter(order, fchp/(fs/2), 'high'); % 'high' specifies a high-pass filter

% Filter the magnitude of acceleration
filteredMag = filtfilt(b1, a1, mag);

% Plot raw vs filtered acceleration
figure;
plot(time, mag); hold on;
plot(time, filteredMag, 'r');
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
legend('Raw', 'Filtered');
title('Acceleration Data');
hold off;

% First integration (from acceleration to velocity)
velocity = cumtrapz(time, filteredMag);

% You might want to apply another high-pass filter to the velocity if there's still a drift
[b2, a2] = butter(order, fchp/(fs/2), 'high');
filteredVelocity = filtfilt(b2, a2, velocity);

% Plot the velocity
figure;
plot(time, filteredVelocity);
xlabel('Time (s)');
ylabel('Velocity (m/s)');
title('Velocity Data');
axis tight;

% Second integration (from velocity to displacement)
displacement = cumtrapz(time, filteredVelocity);

% Again, you might want to apply a high-pass filter to the displacement if needed
[b3, a3] = butter(order, fchp/(fs/2), 'high');
filteredDisplacement = filtfilt(b3, a3, displacement);

% Plot the displacement
figure;
plot(time, filteredDisplacement);
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Displacement Data');
axis tight;

------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

[Stella]

Thank you for your answer!! I also tried your code and I may use it on my other data!! I really appreciate your fast response!!!