Azzera filtri
Azzera filtri

Filter a noisy signal

2 visualizzazioni (ultimi 30 giorni)
Merkhav E
Merkhav E il 9 Mar 2021
Commentato: Mathieu NOE il 12 Mar 2021
Hello, I have data from a pair of strain gages, and is very noisy. How can I filter the data to get a cleaner plot and be able to identify the peaks?
Thank you in advance for your help.
load('Data_8.mat');
R=T.Right;
L=T.Left;
t=T.X_value;
figure
plot(t,R);hold on;plot(t,L);legend('Right','Left');
grid

Risposta accettata

Mathieu NOE
Mathieu NOE il 9 Mar 2021
hello
you can use many ways to lowpass filter your data
existing matlab functions to smooth noisy data is smooth or movmean but I prefer my slidingavg function
here 's the code :
load('Data_8.mat');
R=T.Right;
L=T.Left;
t=T.X_value;
figure(1)
plot(t,R);hold on;plot(t,L);legend('Right','Left');
grid
figure(2)
N = 2500;
Rs = slidingavg(R, N);
Ls = slidingavg(L, N);
plot(t,Rs);hold on;plot(t,Ls);legend('Right','Left');
grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 michele@giugliano.info)
%
% 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
  2 Commenti
Merkhav E
Merkhav E il 9 Mar 2021
Thank you @Mathieu NOE, it works good. How did you select the value for N =2500?
Mathieu NOE
Mathieu NOE il 9 Mar 2021
A compromise between smoothing the noisy segments while preserving the useful peaks and steps edges sharp enough
I noticed that your sampling frequency is pretty high and that I could decimate the signals by a factor as high as 50 without changing too much the shape of the signal
This can be a benefit for speeding up the smoothing process as you reduce by a significant factor the amount of data to process
a code to test the different approaches - including suggestions made by Image Analyst
load('Data_8.mat');
R=T.Right;
L=T.Left;
t=T.X_value;
samples = length(t);
Fs = (samples-1)/(t(samples)-t(1));
figure(1)
plot(t,R);hold on;plot(t,L);legend('Right','Left');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz' ]);
grid
% NB : decim = 1 will do nothing (output = input)
decim = 50;
if decim>1
R = decimate(R,decim);
L = decimate(L,decim);
Fs = Fs/decim;
end
samples = length(R);
t = (0:samples-1)*1/Fs;
figure(2)
plot(t,R);hold on;plot(t,L);legend('Right','Left');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz' ]);
grid on
figure(3)
N = 25;
Rs = slidingavg(R, N);
Ls = slidingavg(L, N);
plot(t,Rs);hold on;plot(t,Ls);legend('Right','Left');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with slidingavg' ]);
grid on
figure(4)
N = 50;
Rs = medfilt1(R, N,'truncate');
Ls = medfilt1(L, N,'truncate');
plot(t,Rs);hold on;plot(t,Ls);legend('Right','Left');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with medfilt1' ]);
grid on
figure(5)
N = 50;
Rs = sgolayfilt(R,3,51);
Ls = sgolayfilt(L,3,51);
plot(t,Rs);hold on;plot(t,Ls);legend('Right','Left');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Smoothed with sgolayfilt' ]);
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 michele@giugliano.info)
%
% 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

Accedi per commentare.

Più risposte (2)

Image Analyst
Image Analyst il 9 Mar 2021
You can use a median filter if you want to smooth the noisy parts but keep the step edges sharp. Use medfilt1() in the Signal Processing Toolbox.
You can also fit it with a sliding polynomial fit with sgolayfilt(). Using that with a polynomial order of 2 or 3 will let you follow and keep true peaks better than a sliding average filter (which is a polynomial order of 1).

Merkhav E
Merkhav E il 10 Mar 2021
Excellent. Thank you again. I can show the max value, but how can I show not just the max, but let´s say the 5 highest peaks or the min values. I tried [pks,lcs,w,p] = findpeaks (...), but not getting the plot.
load('Data_8.mat');
R=T.Right;
L=T.Left;
t=T.X_value;
samples = length(t);
Fs = (samples-1)/(t(samples)-t(1));
[R_max,idx]= max(R);
[L_max,idx]= max(L);
t_m = t(idx);
figure(1)
plot(t,R);hold on;plot(t,L);legend('Right','Left');
plot(t_m, R_max,'^r');hold on; plot (t_m, L_max,'ro')
text(t_m,R(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', R_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
text(t_m,L(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', L_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
title(['Data samples at Fs = ' num2str(round(Fs)) 'Hz' ]);
grid
% NB : decim = 1 will do nothing (output = input)
decim = 50;
if decim > 1
R = decimate (R,decim);
L = decimate (L,decim);
Fs = Fs/decim;
end
samples = length(R);
t = (0:samples - 1)*1/Fs;
[R_max,idx]= max(R);
[L_max,idx]= max(L);
t_m = t(idx);
figure(2)
plot(t,R); hold on; plot (t,L); legend('Right','Left');
plot(t_m, R_max,'^r');hold on; plot (t_m, L_max,'ro')
text(t_m,R(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', R_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
text(t_m,L(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', L_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
title(['Data samples at Fs = ' num2str(round(Fs)) 'Hz']);
grid on
figure(3)
N = 25;
Rs = slidingavg(R,N);
Ls = slidingavg(L,N);
[R_max,idx]= max(Rs);
[L_max,idx]= max(Ls);
t_m = t(idx);
plot(t,Rs); hold on; plot (t,Ls); legend('Right','Left');
plot(t_m, R_max,'^r');hold on; plot (t_m, L_max,'ro')
text(t_m,R(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', R_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
text(t_m,L(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', L_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
title(['Data samples at Fs = ' num2str(round(Fs)) 'Hz / Smoothed with slidingavg' ]);
grid on
figure(4)
N = 50;
Rs = medfilt1(R, N, 'truncate');
Ls = medfilt1(L, N, 'truncate');
[R_max,idx]= max(Rs);
[L_max,idx]= max(Ls);
t_m = t(idx);
plot(t,Rs); hold on; plot (t,Ls); legend('Right','Left');
plot(t_m, R_max,'^r');hold on; plot (t_m, L_max,'ro')
text(t_m,R(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', R_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
text(t_m,L(idx), sprintf('\\leftarrow Max = %.6f\n t = %.2f ', L_max,t_m), 'HorizontalAlignment','left', 'VerticalAlignment','top')
title(['Data samples at Fs = ' num2str(round(Fs)) 'Hz / Smoothed with medfilt1' ]);
grid on
figure(5)
N = 50;
Rs = sgolayfilt(R,3,51);
Ls = sgolayfilt(L,3,51);
[Rs_max,index]= max(Rs);
[Ls_max,index]= max(Ls);
t_max = t(index);
plot(t,Rs); hold on; plot (t,Ls); legend('Right','Left');
plot(t_max, Rs_max,'^r');hold on; plot (t_max, Ls_max,'ro')
text(t_max,Rs(index),sprintf('\\leftarrow Max = %.6f\n t = %.2f ', Rs_max, t_max), 'HorizontalAlignment','left', 'VerticalAlignment','top')
text(t_max,Ls(index),sprintf('\\leftarrow Max = %.6f\n t = %.2f ', Ls_max, t_max), 'HorizontalAlignment','left', 'VerticalAlignment','top')
title(['Data samples at Fs = ' num2str(round(Fs)) 'Hz / Smoothed with sgolayfilt' ]);
grid on
  7 Commenti
Image Analyst
Image Analyst il 12 Mar 2021
In text() you can use the 'VerticalAlignment' option for 'bottom'.
Mathieu NOE
Mathieu NOE il 12 Mar 2021
tx but it the same messy display anyway

Accedi per commentare.

Categorie

Scopri di più su Descriptive Statistics in Help Center e File Exchange

Tag

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by