Remove outliers but be careful with end points

Hi, I have some typcial data like this that I want to remove the outliers (red arrows)

I use the rmoutliers function and then fillmissing to handle these.

        %For last plotted Data
        ax=app.UIAxes3;
        [datax,datay] = getDataFromGraph(app,ax,1);   % my fucntion that gets the last plot x and y values.
     
        A=[datax datay];
        [B,TFrm,TFoutlier] = rmoutliers(A,"movmedian",10);
        x=B(:,1); y=B(:,2);     hold(ax,"on");

The contraint I have is if I remove an outlier from the 1st or last data point (i.e. green arrow above), then it MUST be replaced with e.g the nearest non outlier.

So I thought this would do it:

        F = fillmissing(y,'linear','EndValues','nearest');       %F = fillmissing(y,'movmedian',10);  
        plot(ax,x,F,'.-'); 

However, its ignoring the last point (I dont mind other outliers being ignored, I just need the starting x and finishing x to be the same as the original data

Thanks

Answer 1

[Star Strider]

The rmoutliers function does not know that the end points are outliers because it has nothing after them to compare with beyond that. (Neither do you, actually. They could be valid data.)

You could supply an additional end value (perhaps the mean of the previous values), or just remove them yourself.

EDIT — Corrected typographical errors.

[Jason]

But how would I know if the end point was actually an outlier?

[Star Strider]

That's the point! You can't know, and neither does the function.

There ar many different ways of removing noise. The Savitzky-Golay filter is one.

Another option is a simple lowpass IIR filter —

A1 = readmatrix('myData.csv')

A1 = 100×2

95.0000 0 95.4000 0.0000 95.6000 -0.0000 95.8000 -0.0000 96.0000 0.0000 96.2000 0.0000 96.4000 -0.0001 96.6000 -0.0000 96.8000 -0.0000 97.0000 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = A1(:,1);
s = A1(:,2);

figure
plot(t, s)
grid
title('Original Data')

[FTs1,Fv] = FFT1(s,t);

[vys, vlocs] = findpeaks(-abs(FTs1)*2);



figure
plot(Fv, abs(FTs1)*2)
grid
xline(Fv(vlocs(1)), '--', 'F_{co}')
xlabel('Frequency')
ylabel('Magnitude')

●

Fs = 1/mean(diff(t));
Fco = Fv(vlocs(1));

s_filt = lowpass(s, Fco, Fs, ImpulseResponse='iir');

figure
plot(t, s, DisplayName='Original Data')
hold on
plot(t, s_filt, DisplayName=sprintf('Lowpass Filtered Fco = %.3f',Fco))
hold off
grid
legend(Location='NW')
title('Original & Lowpass-Filtered Data')

function [FTs1,Fv] = FFT1(s,t)
% One-Sided Numerical Fourier Transform
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

You can set the cutoff frequency 'Fco' to be whatever you want. Here, is set it to the first valley in the fft result.

.

[Jason]

Thankyou for this

[Jason]

Is there a reason why this doesn't work:

ax=app.UIAxes3;
[datax,datay] = getDataFromGraph(app,ax,1);
head(datax)
head(datay)
[FTs1,Fv] = FFT1(app,datay,datax);
[vys, vlocs] = findpeaks(-abs(FTs1)*2);
Fs = 1/mean(diff(datax));
Fco = Fv(vlocs(1));
s_filt = lowpass(datay, Fco, Fs, ImpulseResponse='iir');
% figure
% plot(t, s, DisplayName='Original Data')
hold(ax,'on');
plot(ax,datax, s_filt, DisplayName=sprintf('Lowpass Filtered Fco = %.3f',Fco))
hold(ax,'off')

title(ax,'Original & Lowpass-Filtered Data')

head(datax) and head(datay) gives:

           95
         95.4
         95.6
         95.8
           96
         96.2
         96.4
         96.6
         
            0
   1.3378e-05
  -3.2534e-05
  -3.9384e-06
   1.6481e-06
   1.1045e-05
   -6.113e-05
  -5.5223e-06

the error I'm getting is:

  Undefined function 'hann' for input arguments of type 'double'.
Error in HTS_TestSoftware/LowPassButtonPushed (line 13018)
            [FTs1,Fv] = FFT1(app,datay,datax);

[Jason]

Ahh... I see the signal toolbox is required which I don't have!

hann requires Signal Processing Toolbox.

[Star Strider]

As always, my pleasure!

When you mentioned the Savitzky-Golay filter, I thought you had the Signal Processing Toolbox. I was under the impression that it was required for the smoothdata function 'sgolay' option, since I encountered that problem before with someone who could not use the 'sgolay' option because of they did not have the Signal Processing Toolbox. That must have changed since.

You can do the filtering here, including using my 'FFT1' function, save the results (to a .txt or .csv file or something else), and then download that file from MATLAB Drive either directly oir using the websave function:

file = websave('filename.csv', 'https://URL path/filename.csv');

You only have to know the URL of the file.

For example, to get the file you attached —

file = websave('myData.csv', 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1833407/myData.csv')

file = '/users/mss.system.r858k/myData.csv'

Data = readtable(file)

Data = 100x2 table
    Var1       Var2    
    ____    ___________

      95              0
    95.4     1.3378e-05
    95.6    -3.2534e-05
    95.8    -3.9384e-06
      96     1.6481e-06
    96.2     1.1045e-05
    96.4     -6.113e-05
    96.6    -5.5223e-06
    96.8    -3.7543e-05
      97    -4.9144e-05
    97.2    -3.2042e-05
    97.4     5.2269e-05
    97.6    -2.7162e-05
    97.8    -4.7838e-05
      98    -2.4144e-05
    98.2    -4.4071e-05

I could do that directly in my MATLAB installation on this computer as well as here in Answers. The file would be imported and I would have access to it. I could then write it to a local directory on my computer.

.

[Jason]

Ah thats a shame. This is part of an instrument control system where i need almost instant smoothing and I cant be accessing the web to get results

[Star Strider]

If any part of your system is analog, your best option then may be to design a Bessel filter with the appropriate lowpass characteristic and realise it in analog hardware, for example just ahead of the ADC stage. (Bessel filters are phase-neutrral, so there is no phase distortion or phase delay.) ICs for this purpose exist, although I've never used them.

If you have the Control System Toolbox, you might be able to design an appropriate continuous-time filter (using MATLAB Answers) and implement the filter in discrete time with a state-space or transfer-function realisation of it. (Bessel filters lose thier phase-neutral characteristic if implemented as discrete filters. An elliptic filter would be my choice.)

Your best option however is probably to get the Signal Processing Toolbox, although the 'almost real time' constraint may prohibit that working as well as you might need it to work.

Answer 2

[Steven Lord]

Do you want to remove the outliers or do you want to fill them in? If the latter, see the filloutliers function.

[Jason]

Well ideally remove them except any that are end points:

I have this so far,

B = filloutliers(A,'center','movmedian',10);
x=B(:,1); y=B(:,2);     hold(ax,"on");
plot(ax,x,y,'y.-'); 

but its not quite what I want

Ultimately I want a nice smooth curve like this red one that I obtained by performing sGolay smooting on - I just need it to span the whole x-range as the original data

(Data attached)

[Steven Lord]

If you want to smooth your data using Savitzky-Golay, you can use the smoothdata function and specify the "sgolay" smoothing method.

[Jason]

Yes thats I currently do, thanks

degree=app.SGDegreeEditField.Value;        
ys = smoothdata(y,'sgolay','Degree',degree);       %sgolay, Degree default is 2
hold(ax,'on');
plot(ax,x,ys,'y-','LineWidth',0.5); %grid(ax,"on");