How do I take an average around unusual values?

Hi,

I have a matrix with two rows - one that has values of interest and then the other that has corresponding moving standard deviation values. Sometimes there are spikes in the second array (std gets too high) - and there could be several of them. I already know how to identify these spikes - they are 3 standard deviation away from the mean (of moving standard deviation values).

My problem is to remove these spikes once they are identified, by averaging the values in the first row: take a value right before the spike occurred, add the value right after the spike ended (std values are back to being less than 3 standard dev. from the mean), and then divide those by two - and fill the elements corresponding to the spike with those averages.

Could anyone please help? I'm confused with how the indexing would work in this case, since there could be multiple occurrences of an event, and the averages should be taken for each of those events, properly.

Thanks!!

Answer 1

[Jan]

With interp1:

signal  = rand(1000, 1);  % Test data, the first column of your matrix
isSpike = rand(1000, 1) < 0.01;
signal(isSpike) = interp1(signal(~isSpike), find(~isSpike), find(isSpike));

Answer 2

[Image Analyst]

If you have the Image Processing Toolbox, you could use roifill(). Otherwise you could identify the "good" elements, pass them into interp1(), and get estimates for the "bad/missing" areas. If you want a higher order estimate (e.g. quadratic) then you could filter the signal with a Savitzky-Golay filter (sgolay() in the Signal Processing Toolbox) and use the filtered values to replace the values. Something like (untested):

badIndexes = stdDevSignal > someValue;
filteredSignal = sgolay(originalSignal, 2); % or use roifill() or interp1()
goodSignal = originalSignal; % Initialize
goodSignal(badIndexes) = filteredSignal(badIndexes); % Replace bad values.

[Jan]

And if you do not have the Signal-Processing-Toolbox or want more speed: FEX: fSgolayFilt.

Answer 3

[Veronika]

Thanks! I see that the interp1 turns those values corresponding to a spike into NaN values. Is there a way to instead replace them with the local average values?

I've tried so far nanmean and inpaint_nans - they seem to help with the problem, but do not really provide local averages.

[Jan]

INTERP1 does not insert NaNs for Spikes, except they are found on the edges. So please post the code, which creates the NaNs in your case.

[Veronika]

spike = s > mean(s) + 3*std(s);

sig(spike) = interp1(sig(~spike), find(~spike), find(spike));

sig = inpaint_nans(sig);