Best filter for velocity noisy data

Hello,

First of all thank you for your time considering my doubt. I am not a programmer nor have much idea of matlab...

I have a vector which is part of a speed signal with too much noise. I have tried several filtering mehods but the best I can do has an important lag on the data. Is there any possibility to filter the data without that lag?

You can see what I mean on the figure attached.

About the filter I am using, the code is:

    % Filter speed
    fc = 1;
    fs = 100;
    Wn = fc/(fs/2);
    n = 2;
    [b,a] = butter(n, Wn);
    speed = filtfilt(b, a, raw_speed);
    speed = 2*pi*speed;

For me the shape of the curve is ok, regardless the moments when the raw speed is touching the zero line, but not the filter signal.

Find also attached the data.

Thank you.

[Mathieu NOE]

hello

lag is a time shift , but as you use filtfilt (good point ) and not filter, the phase lag of the butter filetr is cancelled.

So there is "no" time lag between the raw and filtered signal

all other "smoothing" techniques available will have more less the same result - see example below.

I thought that your concern would have been more about the filtered signal being too different in the areas where it touches zero.

load('matlab.mat');

y = raw_speed_rads;

Fs = 100;
samples = length(y);
dt = 1/Fs;
t = (0:samples-1)*dt;

% %%%%%%%%%%%%%%%%
figure(1)
N = 25;
ys = slidingavg(y, N);
plot(t,y,t,ys);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with slidingavg' ]);

% %%%%%%%%%%%%%%%%
figure(2)
N = 25;
ys = medfilt1(y, N,'truncate');
plot(t,y,t,ys);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with medfilt1' ]);
grid on

%%%%%%%%%%%%%%%%
figure(3)
ys = sgolayfilt(y,5,51);
plot(t,y,t,ys);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with sgolayfilt' ]);
grid on

%%%%%%%%%%%%%%%%
NN = 4;
Wn = 0.1;
[B,A] = butter(NN,Wn);

figure(4)
ys = filtfilt(B,A,y);
plot(t,y,t,ys);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with butterworth LP' ]);
grid on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = slidingavg(in, N)
%   OUTPUT_ARRAY = SLIDINGAVG(INPUT_ARRAY, N)
%
%  The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
%  the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
%  vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
%  of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
%  *  The input argument INPUT_ARRAY is the numerical data array to be processed.
%  *  The input argument N  is the number of neighboring data points to average over for each point of IN.
%
%  *  The output argument OUTPUT_ARRAY is the output data array.
%
%     © 2002 - Michele Giugliano, PhD and Maura Arsiero
%     (Bern, Friday July 5th, 2002 - 21:10)
%    (http://www.giugliano.info) (bug-reports to [e-mail])
%
% Two simple examples with second- and third-order filters are
%  slidingavg([4 3 5 2 8 9 1],2)   
%  ans =
%   3.5000  4.0000  3.3333  5.0000  6.3333  6.0000  5.0000
%
%  slidingavg([4 3 5 2 8 9 1],3)
%  ans =
%   3.5000  4.0000  3.3333  5.0000  6.3333  6.0000  5.0000
%

if (isempty(in)) | (N<=0)                                              % If the input array is empty or N is non-positive,
 disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));     % an error is reported to the standard output and the
 return;                                                               % execution of the routine is stopped.
end % if

if (N==1)                                                              % If the number of neighbouring points over which the sliding 
 out = in;                                                             % average will be performed is '1', then no average actually occur and
 return;                                                               % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if                                                               % is stopped.

nx   = length(in);             % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.

if (N>=(2*(nx-1)))                                                     % If the number of neighbouring points over which the sliding 
 out = mean(in)*ones(size(in));                                        % average will be performed is large enough, then the average actually covers all the points
 return;                                                               % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if                                                               % The execution of the routine is stopped.



out = zeros(size(in));         % In all the other situations, the initialization of the output data structure is performed.

if rem(N,2)~=1                 % When N is even, then we proceed in taking the half of it:
 m = N/2;                      % m = N     /  2.
else                           % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
 m = (N-1)/2;                  % m = (N-1) /  2.
end % if

for i=1:nx,                                                 % For each element (i-th) contained in the input numerical array, a check must be performed:
 if ((i-m) < 1) & ((i+m) <= nx)                             % If not enough points are available on the left of the i-th element..
  out(i) = mean(in(1:i+m));                                 % then we proceed to evaluate the mean from the first element to the (i + m)-th.
 elseif ((i-m) >= 1) & ((i+m) <= nx)                        % If enough points are available on the left and on the right of the i-th element..
  out(i) = mean(in(i-m:i+m));                               % then we proceed to evaluate the mean on 2*m elements centered on the i-th position.
elseif ((i-m) >= 1) & ((i+m) > nx)                          % If not enough points are available on the rigth of the i-th element..
  out(i) = mean(in(i-m:nx));                                % then we proceed to evaluate the mean from the element (i - m)-th to the last one.
 elseif ((i-m) < 1) & ((i+m) > nx)                          % If not enough points are available on the left and on the rigth of the i-th element..
  out(i) = mean(in(1:nx));                                  % then we proceed to evaluate the mean from the first element to the last.
 end % if
end % for i
end

[Alejandro Muñoz López]

Hello, thank you for the reply. You are right, I am sorry bout using the term "lag". My concerns relate to areas arround 0, I need that the filtered signal will include those data comming from 0, because when I do the average of each curve (each curve means a repetition of a strenght exercise, 5 curves equals to 5 reps here) data is not correct.

Any solution to that? Thank you

Answer 1

[Mathieu NOE]

Suggestion of the day

there are certainly other ways to implement a filter that action is function of signal amplitude (or deviation or noise or ...)

I tried first with proportionnal gain , then with non linear action, ended up with simple 3 states filter

clc
close all
clearvars

load('matlab.mat');

y = raw_speed_rads;

Fs = 100;
samples = length(y);
dt = 1/Fs;
t = (0:samples-1)*dt;


%%%%%%%%%%%%%%%%
% smart filter V2 : simple first order low pass filter (2 stages in series), 
% the filtering effect is increased with signal amplitude in 3 discrete steps
% like filtfilt, the filtering is applied once in forward time and back in
% reversed time to cancel phase lag

threshold1 = 0.05;  smoothing_factor1 = 0.65;
threshold2 = 0.5;   smoothing_factor2 = 0.85;

% forward filering (2 stage in series of first order filter)
y_norm = abs(y)./max(y); % signal normalization (for flter parameters usage)
alpha = zeros(size(y));
alpha(y_norm>threshold1) = smoothing_factor1; 
alpha(y_norm>threshold2) = smoothing_factor2; 

ys = zeros(size(y));
ys2 = zeros(size(y));
ys(1) = y(1);
ys2(1) = ys(1);

for ci = 2:samples
    ys(ci) = alpha(ci).*ys(ci-1)  + (1-alpha(ci)).*y(ci);
    ys2(ci) = alpha(ci).*ys2(ci-1)  + (1-alpha(ci)).*ys(ci);
 end
clear ys

% backward filtering (2 stage in series of first order filter)
y_norm = abs(ys2)./max(ys2); % signal normalization (for flter parameters usage)
alpha = zeros(size(y));
alpha(y_norm>threshold1) = smoothing_factor1; 
alpha(y_norm>threshold2) = smoothing_factor2; 

yss = zeros(size(y));
yss(1) = ys2(samples);
yss2(1) = yss(1);

for ci = 2:samples
    yss(ci) = alpha(samples-ci+1).*yss(ci-1)  + (1-alpha(samples-ci+1)).*ys2(samples-ci+1);
    yss2(ci) = alpha(samples-ci+1).*yss2(ci-1)  + (1-alpha(samples-ci+1)).*yss(ci);
end
clear yss
yss2 = yss2(end:-1:1);

figure(5)
% plot(t,y,t,ys2,t,yss2);legend('Raw','Smoothed fw','Smoothed fw + bckward');
subplot(2,1,1),plot(t,y,t,yss2);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with home made filter' ]);
grid on
subplot(2,1,2),plot(t,alpha);legend('alpha');
grid on

Answer 2

[Alejandro Muñoz López]

Thank you for the custom filter, it is awsome. However, I still need the signal more filtered, like butter with n=2 and fc=1 dig, because I have to compute acceleration from speed (accel = Fsamp*[0 diff(speed)];) and as it is right now the acceleration signal is crazy

Thank you again

[Mathieu NOE]

hello again - more stuff regarding non linear filtering

check this !

probably the best I could do . see also the attaced function for a less noisy derivative

clc
close all
clearvars

load('matlab.mat');
y = raw_speed_rads;
Fs = 100;
samples = length(y);
dt = 1/Fs;
t = (0:samples-1)*dt;

% tre trick here is 
% 1/ to convert to log scale (like mu law encoder) and resample the data at higher frequency 
% 2/ do the non linear filtering
% 3/ and last come back to original linear scale

ylog= log10(y+1e-20);

%% oversampling
factor = 10;
samples = samples*factor;
Fs = Fs*factor;
tt = linspace(min(t), max(t),samples);
ylog=interp1(t,ylog,tt);

%% alpha law
alpha=ylog+30;  % a little bit of shift for dynamic compression (the higher the shift the closer to 1 remains alpha)
alpha=movmean(alpha,5*factor); % a little bit of smoothing to keep low alpha valleys wide enough
alpha=0.995*alpha./max(alpha);
alpha(alpha<0) = 0;


%% %%%%%%%%% main loop %%%%%
ys = zeros(size(y));
ys(1) = y(1);

for ci = 2:samples
    ys(ci) = alpha(ci).*ys(ci-1)  + (1-alpha(ci)).*ylog(ci);
 end

% backward filtering (2 stage in series of first order filter)
yss = zeros(size(y));
yss(1) = ys(samples);

for ci = 2:samples
    yss(ci) = alpha(samples-ci+1).*yss(ci-1)  + (1-alpha(samples-ci+1)).*ys(samples-ci+1);
end
yss = yss(end:-1:1);

%% back conversion
yss = 10.^yss; % converting back from log to linear scale
yss = downsample(yss,factor); % converting back to original sample rate
alpha = downsample(alpha,factor); % converting back to original sample rate

figure(5)
subplot(2,1,1),plot(t,y,t,yss);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with home made filter' ]);
grid on
subplot(2,1,2),plot(t,alpha);legend('alpha');
grid on

%% firt and second derivative
[dy, ddy] = firstsecondderivatives(t,yss);

dys = sgolayfilt(dy,1,21); % more smoothing ?

figure(6)
plot(t,dy,t,dys);
title('Acceleration');legend('Raw','Smoothed');
grid on

[Alejandro Muñoz López]

Thank you again! This could suit for me! However, when I am using another file (attached), it does not work when calculating "ylog= log10(y+1e-20);" because it is creating like a double vector.

Any idea? Thank you in advance again!

[Mathieu NOE]

hello

I found out that this new data had some slightly negative values

that causes the log10 to generate complex and not only real values data

it is now fixed in the code below :

clc
close all
clearvars

% load('matlab.mat');
load('Set3.mat');

y = raw_speed_rads;

Fs = 100;
samples = length(y);
dt = 1/Fs;
t = (0:samples-1)*dt;

% tre trick here is 
% 1/ to convert to log scale (like mu law encoder) and resample the data at higher frequency 
% 2/ do the non linear filtering
% 3/ come back to original linear scale
% 4/ downsample to come back to original sampling frequency

% compensate for slightly negative values (like y(y<0) =   1.0e-13 *   (-0.3906   -0.3906   -0.3906   -0.3906   -0.3906   -0.3906   -0.3906   -0.0663)
y = y-min(y);

ylog= log10(y+1e-20);

%% oversampling
factor = 10;
samples = samples*factor;
Fs = Fs*factor;
tt = linspace(min(t), max(t),samples);
ylog=interp1(t,ylog,tt);

%% alpha law
alpha=ylog+30;  % a little bit of shift for dynamic compression (the higher the shift the closer to 1 remains alpha)
alpha=movmean(alpha,5*factor); % a little bit of smoothing to keep low alpha valleys wide enough
alpha=0.995*alpha./max(alpha);
alpha(alpha<0) = 0;


%% %%%%%%%%% main loop %%%%%
ys = zeros(size(y));
ys(1) = y(1);

for ci = 2:samples
    ys(ci) = alpha(ci).*ys(ci-1)  + (1-alpha(ci)).*ylog(ci);
 end

% backward filtering (2 stage in series of first order filter)
yss = zeros(size(y));
yss(1) = ys(samples);

for ci = 2:samples
    yss(ci) = alpha(samples-ci+1).*yss(ci-1)  + (1-alpha(samples-ci+1)).*ys(samples-ci+1);
end
yss = yss(end:-1:1);

%% back conversion
yss = 10.^yss; % converting back from log to linear scale
yss = downsample(yss,factor); % converting back to original sample rate
alpha = downsample(alpha,factor); % converting back to original sample rate

figure(5)
subplot(2,1,1),plot(t,y,t,yss);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with home made filter' ]);
grid on
subplot(2,1,2),plot(t,alpha);legend('alpha');
grid on

%% firt and second derivative
[dy, ddy] = firstsecondderivatives(t,yss);

dys = sgolayfilt(dy,1,21); % more smoothing ?

figure(6)
plot(t,dy,t,dys);
title('Acceleration');legend('Raw','Smoothed');
grid on

[Alejandro Muñoz López]

Now it is perfect! thank you!

[Mathieu NOE]

you're welcome