matlab - Detecting & Removing steps in signal

Given the signal shown below, what is the best way to remove the steps and local maximas it contains. The signal contains some steps which can last up to 100 Samples before they return to about the same value as before the step (Marked with red circles). There are also some peaks which last only for one sample which should also be removed.

First Dataset: On this set, the median Filter works as expected (not shown here).

On the following Dataset, the Medianfilter with the same settings fails. To many samples are replaced by the median value.

What I've tried already:

• Work with a moving average. Advantage: Small peaks have less influence and are nicely removed. Disadvantage: Long lasting steps remain in the signal.


• Get the first difference of the signal (diff(X)) and correct the value in X if a big positive change is followed by a negative change immediately. Advantages and disadvatages are the same as above.

What are some good methods to remove/correct these unwanted values?

Answer

The problem with the moving average is that the average is not robust to the presence of the outliers - so you would need a very large window size to "dilute" the outliers.

Try a non-linear filter instead, like a median filter: Apply a median filter on your signal - you would need a window size of at least 300 samples. Compute the difference between the original signal and the median-filtered version. If the difference is above a threshold, replace the signal by the median-filtered version.

---

Here is some scilab code that attempts to implement this suggestion. The results are plotted here; it seems to work nicely on the fabricated data.

function sm = smooth(x,len)
    sm = filter(ones(1,len)/len,1,x);
endfunction

N = 6000;
x = 0.1*rand(1,N,'normal');
y = cumsum(x);

y = smooth(y,100);

clf
subplot(211);
plot(y)


Njumps = 20;

jump_indices = round(rand(1,Njumps)*N);


jump_length = 0;

y2 = y;

for idx = jump_indices,
    y2(min(N,idx:(idx+jump_length))) = y(min(N,idx:(idx+jump_length))) + 1;
    //plot(min(N,idx:(idx+jump_length)),y2(min(N,idx:(idx+jump_length))),'r')
    jump_length = jump_length + 1;
end


plot(y2,'g');

filter_length = 10;
y3 = y;
for k = 1:N,
    y3(k) = median(y3(max(1,min(N,(k-filter_length/2):(k+filter_length/2)))));
end

plot(y3,'k')


subplot(212);
plot(y-y2);
plot(y-y3,'r');