xdata = (0:0.1:2*pi)';
y0 = sin(xdata);

Add noise to the signal with nonconstant variance.

% Response-dependent Gaussian noise
gnoise = y0.*randn(size(y0));

% Salt-and-pepper noise
spnoise = zeros(size(y0));
p = randperm(length(y0));
sppoints = p(1:round(length(p)/5));
spnoise(sppoints) = 5*sign(y0(sppoints));

ydata = y0 + gnoise + spnoise;

Fit the noisy data with a baseline sinusoidal model, and specify 3 output arguments to get fitting information including residuals.

f = fittype('a*sin(b*x)');
[fit1,gof,fitinfo] = fit(xdata,ydata,f,'StartPoint',[1 1]);

Examine the information in the fitinfo structure.

fitinfo

fitinfo = struct with fields:
           numobs: 63
         numparam: 2
        residuals: [63x1 double]
         Jacobian: [63x2 double]
         exitflag: 3
    firstorderopt: 0.0883
       iterations: 5
        funcCount: 18
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 4.1539e-04
          message: 'Success, but fitting stopped because change in residuals less than tolerance (TolFun).'

Get the residuals from the fitinfo structure.

residuals = fitinfo.residuals;

Identify "outliers" as points at an arbitrary distance greater than 1.5 standard deviations from the baseline model, and refit the data with the outliers excluded.

I = abs( residuals) > 1.5 * std( residuals );
outliers = excludedata(xdata,ydata,'indices',I);

fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],...
           'Exclude',outliers);

Compare the effect of excluding the outliers with the effect of giving them lower bisquare weight in a robust fit.

fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');

Plot the data, the outliers, and the results of the fits. Specify an informative legend.

plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*')
hold on
plot(fit2,'c--')
plot(fit3,'b:')
xlim([0 2*pi])
legend( 'Data', 'Data excluded from second fit', 'Original fit',...
    'Fit with points excluded', 'Robust fit' )
hold off