Noisy plot after deviation (no sensor)

Question

Bruce Rogers 2021-6-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/851865-noisy-plot-after-deviation-no-sensor

评论： Mathieu NOE 2021-6-15

采纳的回答： Mathieu NOE

z_irl_data.txt

在 MATLAB Online 中打开

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Radhey Shyam Meena 2021-6-9

编辑：Radhey Shyam Meena 2021-6-9

在 MATLAB Online 中打开

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.

Bruce Rogers 2021-6-10

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 2021-6-10

@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 2021-6-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/851865-noisy-plot-after-deviation-no-sensor#answer_721605

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2021-6-14

在 MATLAB Online 中打开

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
% brueckner@sbrs.net
% 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 2021-6-15

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 2021-6-15

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 2021-6-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/851865-noisy-plot-after-deviation-no-sensor#answer_721530

在 MATLAB Online 中打开

@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');

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Bruce Rogers 2021-6-10

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.

请先登录，再进行评论。