Noisy plot after deviation (no sensor)

Hello,

i'm trying to get the acceleration from position and time data. (I dont get the data from an accelerometer, this would explain the noise! )The plot seems to be really noisy, is there a good explanation? Shouldn't I be able to see a perfect sine wave? I'm thankful for every idea, thanks!

Here my code and a picture of the result

h = 1:height(z_irl);

for i = 1:height(z_irl)

    t(i) = 0.004 * h(i);

end

t = t';

dt = diff(t); 
dz_irl = diff(z_irl);
%velocity

vel_irl = dz_irl./dt(1:end);

dvel_irl = diff(vel_irl);
%acceleration
acc_irl = dvel_irl./dt(2:end);

acc_irl(numel(t)) = 0;

figure(4)
plot(t,acc_irl,'-b')
xlabel('time')
ylabel('acceleration')
hold off

[Scott MacKenzie]

Your code doesn't execute --> Unrecognized function or variable 't'.

Perhaps update your question, loading t with the required data.

[Bruce Rogers]

Hello Scott,

I adjusted the question and the code. t is just the time in a distance of 0.004s. I'll upload the z_irl data in a minute

[Radhey Shyam Meena]

Hi Bruce,

I understand that you are trying to plot acceleration over time, after computing the acceleration from the position values as specified in the z_irl_data.txt file.

Assuming, that you accept the values computed for acceleration since you did not raise a concern for their accuracy or authenticity, and your query is more about an explanation of the plot you are seeing.

Since the data is quite large (18k+), and it is squeezed along the time axis, you might not be able to visualize it correctly if you fit the entire data in one figure, and hence making it look like noisy.

One way is to look at only the values of first 50 points, to visualize the data correctly. That can be done as

plot(t(1:50),acc_irl(1:50),'-b')

And verify if this is the visualization you were looking for.

this is how it looks for me:

If you are concerned about the accuracy of the acceleration data OR you want to see a more smooth curve, probably by smoothening the values out on multiple data points (like moving average) OR Otherwise, you can update the question accordingly.

I hope this helps.

Regards.

[Scott MacKenzie]

OK, I played with your code a bit. I don't see any issues. Perhaps the noise in the acceleration data is a natural follow-on to the noise in the velocity data. Below is what I put together. I'm not sure if this is of any help, but here you go, anyway. Good luck.

z_irl = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/647435/z_irl_data.txt');

t = 0.004 * (1:length(z_irl));
t = t';
dt = [0; diff(t)]; 

tiledlayout(3,1);

% position
nexttile;
plot(t,z_irl);
xlabel('time');
ylabel('position');

% velocity
dz_irl = [0; diff(z_irl)];
vel_irl = dz_irl./dt;

nexttile;
plot(t, vel_irl);
xlabel('time');
ylabel('velocity');

% acceleration
dvel_irl = [0; diff(vel_irl)];
acc_irl = dvel_irl./dt;

nexttile;
plot(t, acc_irl, '-b');
xlabel('time');
ylabel('acceleration');

[Bruce Rogers]

Hello @Radhey Shyam Meena,

thanks for your help, but unfortunately the acceleration should look like a sine wave, because the position data shows a negative sine wave. Even if I would only look at 50 points, in total a sine curve can not be detected. Is there some kind of filter that smoothens my curve or can interpolate and calculate the noise?

But many thanks for your help!

[Bruce Rogers]

Hello @Scott MacKenzie,

thanks for proving my code. It looks really strange to me, it's really not the result I expected. I hope there is a way to remove the noise in the velocity data, maybe a filter or something like that.

Thank you for your help!

[Scott MacKenzie]

@Bruce Rogers I'm going to post a solution in a minute that shows the velocity and acceleration as a sine wave.

Answer 1

[Mathieu NOE]

hello

my suggestion, with some smoothing at all stages :

z_irl = readmatrix('z_irl_data.txt');
z_irls = smoothdata(z_irl,'gaussian',100);

dt = 0.004;
t = dt * (1:length(z_irl));
t = t';

tiledlayout(3,1);
% position
nexttile;
plot(t,z_irl, t, z_irls);
xlabel('time');
ylabel('position');

% velocity
vel_irl = gradient(z_irls,dt);
vel_irls = smoothdata(vel_irl,'gaussian',100);
nexttile;
plot(t, vel_irl, t, vel_irls);
xlabel('time');
ylabel('velocity');

% acceleration
acc_irl = gradient(vel_irls,dt);
acc_irls = smoothdata(acc_irl,'gaussian',100);

nexttile;
plot(t, acc_irl, '-b', t, acc_irls);
xlabel('time');
ylabel('acceleration');

I don't think you can really get a sinusoidal acceleration... seems more to be trapezoidal somehow

[Bruce Rogers]

Hello Mathieu,

thanks for your answer. I can't explain, why the acceleration isn't just a sine wave, the second derivate of the position of a negative sine wave should be at least the sine wave. But many thanks, it works for me.

[Mathieu NOE]

hello Bruce

I tried another approach, by using a low order Butterworth bandpass filter (if you increase the order the transients last longer).

First I needed to know the period / frequency of your signal, so that 's where I used the zero crossing funtion crossing_V7.m. This way I can center the bandpass cut off frequencies around the signal frequency;

Then , if you play with the lower and upper cut off frequencies of the filter, you have either a more triangular shape acceleration (if you make the frequency bandwith large so that the upper harmonics are not too much filtered). If you narrow the filter , you can get a sinusoidal acceleration plot, and that is simply because your filter has removed the upper harmonics;

% parameter for "triangular" output (more or less the same output as initial code with smoothdata)

%% bandpass filter section %%%%%%
f_low = 0.5*fc;
f_high = 20*fc;
N = 2;

% parameter for "sinus" output

%% bandpass filter section %%%%%%
f_low = 0.5*fc;
f_high = 2*fc;
N = 2;

% this is now the new code :

z_irl = readmatrix('z_irl_data.txt');
dt = 0.004;
Fs = 1/dt;
t = dt * (1:length(z_irl));
t = t';

threshold = mean(z_irl);
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(z_irl,t,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

period = mean(diff(t0_pos))
fc = 1/period;

% z_irls = smoothdata(z_irl,'gaussian',100);
%% bandpass filter section %%%%%%
f_low = 0.5*fc;
f_high = 20*fc;
N = 2;
[b,a] = butter(N,2/Fs*[f_low f_high]);
% now let's filter the signal
z_irls = filtfilt(b,a,z_irl);
z_irls = z_irls+mean(z_irl);

tiledlayout(3,1);
% position
nexttile;
plot(t,z_irl, t, z_irls);
xlabel('time');
ylabel('position');
% velocity
vel_irl = gradient(z_irls,dt);
%vel_irls = smoothdata(vel_irl,'gaussian',100);
vel_irls = filtfilt(b,a,vel_irl);

nexttile;
plot(t, vel_irl, t, vel_irls);
xlabel('time');
ylabel('velocity');
% acceleration
acc_irl = gradient(vel_irls,dt);
% acc_irls = smoothdata(acc_irl,'gaussian',100);
acc_irls = filtfilt(b,a,acc_irl);

nexttile;
plot(t, acc_irl, '-b', t, acc_irls);
xlabel('time');
ylabel('acceleration');

% the zero crossing funtion crossing_V7.m is as follows :

function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% M Noe
% added positive or negative slope condition

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            slope_sign(ii) = sign(slope);
            t0(ii) = t0(ii) - S(ind(ii)) * slope;
            s0(ii) = level;
        end
    end
end

% extract the positive slope crossing points 
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);

% extract the negative slope crossing points 
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);

end

% % Addition:
% % Some people like to get the data points closest to the zero crossing,
% % so we return these as well
% [CC,II] = min(abs([S(ind-1) ; S(ind) ; S(ind+1)]),[],1); 
% ind2 = ind + (II-2); %update indices 
% 
% t0close = t(ind2);
% s0close = S(ind2);

[Bruce Rogers]

Hello Mathieu, this look great thanks for your help! But I think I have to work with the triangular shaped acceleration, because otherwise the acceleration and velocity are too low as they should be in theory. (v = -25 to 25 mm/s and a = -12.5 to 12.5 mm/s^2)

[Mathieu NOE]

Hello Bruce

ok , I just wanted to show you the options , of course you decide which one fits better your needs

have a good day

Answer 2

[Scott MacKenzie]

@Bruce Rogers I added some filtering using MATLAB's smoothdata function. The result is a pretty good sine wave both for velocity and acceleration. smoothdata "returns a moving average of the elements of a vector using a fixed window length that is determined heuristically". The windows length is not specifiied, but you can play with this using parameters described in the documentation. The result is below. Note that the new waveforms experience both phase shift and attenuation as a result of the filtering. The acceleration waveform is very attenuated, indicating noise as large spikes in the raw data. Anyway, hope this helps. Good luck.

z_irl = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/647435/z_irl_data.txt');

t = 0.004 * (1:length(z_irl));
t = t';
dt = [0; diff(t)]; 

tiledlayout(3,1);

% position
nexttile;
plot(t,z_irl);
xlabel('time');
ylabel('position');

% velocity
dz_irl = [0; diff(z_irl)];
vel_irl = dz_irl./dt;
vel_irl = smoothdata(vel_irl); % filter using moving average

nexttile;
plot(t, vel_irl);
xlabel('time');
ylabel('velocity');

% acceleration
dvel_irl = [0; diff(vel_irl)];
acc_irl = dvel_irl./dt;
acc_irl = smoothdata(acc_irl); % filter using moving average

nexttile;
plot(t, acc_irl, '-b');
xlabel('time');
ylabel('acceleration');

[Bruce Rogers]

This is really nice, I was also reading about the smoothdata function, just when you posted your answer. It's great, but it fakes the result a lot. Like the velocity should go from -25 to 25 and the acceleration from -12.5 to 12.5.

I thought about splitting the velocity vector in smaller vectors of 5 or 10 elements and get the average, so I can get a good signal, but also reduce the noise.