How can I smooth a signal while preserving peaks in MATLAB (R2013a)?

I have a noisy signal with sharp peaks. I would like to smoothe the data while preserving the peaks. How can I do that?

Answer 1

[MathWorks Support Team]

In order to smoothe noisy data while preserving the peaks, you can use the function MSSGOLAY, which smoothes a signal with peaks using least-squares polynomial or MSLOWESS, which smoothes a signal with peaks using nonparametric method.

More information on either of these algorithms can be found in the Documentation:

>> doc mssgolay

>> doc mslowess

[Image Analyst]

The original algorithm by Savitzky and Golay assumes the input vector, X, has uniformly spaced separation units, while mssgolay() also allows one that is not uniformly spaced. mssgolay() is in the Bioinformatics Toolbox . sgolay() is in the Signal Processing Toolbox.

Answer 2

[Image Analyst]

Alternatively you could

1. Find the peaks (their indexes and values) with [peakValues, indexes] = findpeaks(signal).
2. Generate an additional smoothed signal somehow, for example with movmean() or sgolayfilt(). smoothedSignal = movmean(signal, 21)
3. Assign the values of the original signal to the corresponding locations in the smoothed signal to replace the smoothed values with the original peak values.

% Create sample signal with peaks
numPoints = 200;
t = linspace(0, 2*pi, numPoints);
signal = cos(t) + 1;
% Add a small amount of noise to it - small enough that we can smooth it away.
signal = signal + 0.5 * rand(1, length(t));
% Put random peaks on it.
indexes1 = randperm(length(t), 8); % 8 randomly located peaks.
signal(indexes1) = signal(indexes1) + 2;
% Plot it.
subplot(3, 1, 1);
plot(t, signal, 'b.-', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
xlabel('time');
ylabel('signal')
% Find the peaks
[peakValues, indexes2] = findpeaks(signal, 'Threshold', 1);
% Plot peaks over original signal.
hold on;
plot(t(indexes2), peakValues, 'rv', 'LineWidth', 2, 'MarkerSize', 15);
title('Original Signal with Peaks Identified');
% Smooth the signal.  (With a little more effort you could smooth only the non-peak values so the peaks don't influence the smoothed result.)
nonPeakindexes = setdiff(1:length(t), indexes2);
smoothedSignal = movmean(signal, 21);
% Plot it.
subplot(3, 1, 2);
plot(t, smoothedSignal, 'b.-', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
xlabel('time');
ylabel('signal')
title('Smoothed Signal');

% Replace the original peak values at their original location and height.
finalSignal = smoothedSignal; % Initialize
finalSignal(indexes2) = signal(indexes2); % Replace smoothed value with original value.
% Plot it.
subplot(3, 1, 3);
plot(t, finalSignal, 'b.-', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
xlabel('time');
ylabel('signal')
title('Smoothed Signal With Peaks Back In');