How to omit narrow variances in data

Dear MATLAB Community,

I have some data, obtained experimentally, that I'd like to refine. In the attached plot, you'll notice that there are some very narrow-band, noisy oscillations between +/- y values. For example, the region x = 920 to 980, the region x = 560 to 640, and the region x = 10 to 60 contain lots of noisy oscillations while the region x = 340 to 450 reveals a smooth, clear, non-noisy curve.

What MATLAB code could I impliment to remove these noisy oscillations? The "smoothdata()" function thus far has been unsucessful, so I'm thinking some kind of an "if" statement would be helpful here. I'm imagining some code that measures the "width" of each oscillation and essentially ignores narrow oscilations.

PS: For reference, the maximum and minimum y values of this plot are +/- radians, and the x axis is in hertz (frequency).

Thank you in advance for your help!

[Walter Roberson]

movmedian maybe ?

[Mathieu NOE]

hello

could you share your data ?

your data contains only +/- pi/2 values or any values between these two extrema ?

what logic would you implement in case of "dense" oscillations : do we put these data to zero ?

Answer 1

[Mathieu NOE]

seems to me that smoothdata with movmedian method is appropriate , which in fact is also what @Walter Roberson also suggested

% dummy data
n = 1000;
x = 10*(0:n-1)/n;
y = pi/2*square(x);

figure(1)
plot(x,y,'k*'); hold on

% create some "bursts" 
a = 100;
l = 10;
ind = randi([a,a+3*l],l,1);
y(ind) = pi/2*(-1).^(1:l);

% create more "bursts" 
a = 400;
l = 40;
ind = randi([a,a+3*l],l,1);
y(ind) = pi/2*(-1).^(1:l);

ys = smoothdata(y,'movmedian',25);

plot(x,y,'b',x,ys,'r');
legend('clean data','corrupted data','filtered data');

[Mathieu NOE]

FYI, this is another approach (taht seems to me quite overcomplicated, but maybe this could be of some help anyway)

this will simply remove the sections of frequent zigzags with NaNs (now you have to decide if this is what you want or not)

% dummy data
n = 1000;
x = 20*(0:n-1)/n;
y = pi/2*square(x);

% create some "bursts" 
a = 100;
l = 10;
ind = randi([a,a+3*l],l,1);
y(ind) = pi/2*(-1).^(1:l);

a = 400;
l = 40;
ind = randi([a,a+3*l],l,1);
y(ind) = pi/2*(-1).^(1:l);

% signal rate of change
threshold = 0; % 

[t0_pos1,t0_neg1] = find_zc(x,y,threshold);
pretrigger = 0.25; % second
posttrigger = 0.25; % second
freqP = interp1(t0_pos1+pretrigger,1./gradient(t0_pos1),x,'linear',0);
freqN = interp1(t0_neg1-posttrigger,1./gradient(t0_neg1),x,'linear',0);

figure(1)
subplot(2,1,1),
plot(x,y,'b.-',t0_pos1,threshold*ones(size(t0_pos1)),'*r',t0_neg1,threshold*ones(size(t0_neg1)),'*g','linewidth',2,'markersize',12);grid on
xlim([x(1) x(end)]);
title('Signal plot');
legend('signal','PS crossing points','NS crossing points');
xlabel('time (s)')
ylabel('Amplitude')
subplot(2,1,2),plot(x,freqP,'b.-',x,freqN,'r.-','linewidth',2,'markersize',12);grid on
xlim([x(1) x(end)]);
title('Signal rate (frequency)');
legend('PS crossing points rate','NS crossing points rate');
xlabel('time (s)')
ylabel('Hz')

% remove data
f_lim = 1; % signal with rapid changes (above this frequency) will be replaced by NaNs
indP = freqP>f_lim;
indN = freqN>f_lim;
yy = y;
yy(indP) = NaN;
yy(indN) = NaN;

figure(2)
plot(x,y,'b.-',x,yy,'r.-');grid on
title('Signal plot');
xlabel('time (s)')
ylabel('Amplitude')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Mathieu NOE]

hello @Adam Harris

problem solved ?