How to extract and smooth/filter the different peaks of a noisy signal with a decreasing and almost oscillating trend

Hi everybody,

I have a noisy signal with a decreasing and almost oscillating trend and I would like to extract and save all the different peaks of this signal in different vectors. I tried several codes I found on Matlab Central or Stack Overflow, but unfortunately a lot of them require the use of Toolboxes I do not have (Signal Processing, Symbolic, Image Processing...)... I use Matlab 2015b with Simulink and the Control System Toolbox.

My signal looks like this :

And I would like to cut it as shown with the red lines.

First I try to smooth my curve in order to be able to find the different peaks and then divide and extract them. Here ( http://stackoverflow.com/questions/1773542/matlab-filter-noisy-ekg-signal ) is a code I found that sum up different mathematical approaches to smooth a signal. Only the moving average approach works without any Toolbox.

%# load the signal

cs = load('Signal.mat'); % Load the .mat file containing the signal with noise to analyse
S1  = cs.num(:,1);
S2  = cs.num(:,2);

x = S2;
% x = x + randn(1,length(x)).*0.18;

%# plot noisy signal
figure
subplot(911), plot(x), set(gca, 'YLim', [-10 10], 'xtick',[])
title('noisy')

% %# sgolay filter (Signal Processing Toolbox)
% frame = 15;
% degree = 0;
% y = sgolayfilt(x, degree, frame);
% subplot(912), plot(y), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('sgolayfilt')

% %# smooth (Curve Fitting Toolbox)
% window = 30;
% %#y = smooth(x, window, 'moving');
% %#y = smooth(x, window/length(x), 'sgolay', 2);
% y = smooth(x, window/length(x), 'rloess');
% subplot(913), plot(y), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('smooth')

%# moving average filter
window = 15;
h = ones(window,1)/window;
y = filter(h, 1, x);
subplot(914), plot(y), set(gca, 'YLim', [-10 10], 'xtick',[])
title('moving average')

% %# moving weighted window (Signal Processing Toolbox)
% window = 7;
% h = gausswin(2*window+1)./window;
% y = zeros(size(x));
% for i=1:length(x)
%     for j=-window:window;
%         if j>-i && j<(length(x)-i+1) 
%             %#y(i) = y(i) + x(i+j) * (1-(j/window)^2)/window;
%             y(i) = y(i) + x(i+j) * h(j+window+1);
%         end
%     end
% end
% subplot(915), plot( y ), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('weighted window')

% %# gaussian (Statistics and Machine Learning Toolbox)
% window = 7;
% h = normpdf( -window:window, 0, fix((2*window+1)/6) );
% y = filter(h, 1, x);
% subplot(916), plot( y ), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('gaussian')

% %# median filter (Signal Processing Toolbox)
% window = 15;
% y = medfilt1(x, window);
% subplot(917), plot(y), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('median')

% %# filter (Signal Processing Toolbox)
% order = 15;
% h = fir1(order, 0.1, rectwin(order+1));
% y = filter(h, 1, x);
% subplot(918), plot( y ), set(gca, 'YLim', [-10 10], 'xtick',[])
% title('fir1')

% %# lowpass Butterworth filter (DSP System Toolbox - Signal Processing
% Toolbox)
% fNorm = 25 / (Fs/2);               %# normalized cutoff frequency
% [b,a] = butter(10, fNorm, 'low');  %# 10th order filter
% y = filtfilt(b, a, x);
% subplot(919), plot(y), set(gca, 'YLim', [-10 10])
% title('butterworth')

The result I get using the moving average method is still noisy and I do not how to cut it in several pieces.

Can anyone advise me regarding this ?

Cheers

Answer 1

[Star Strider]

There are online discrete filter design resources, the best (IMO) being the University of York Butterworth / Bessel / Chebyshev Filters. You can then use the filter function to do the filtering. It has phase distortion (the Signal Processing Toolbox filtfilt function does not), but it will likely work well enough.

Only you know the importance of relevant parts of your data — I don’t — so consider that with respect to the rest of this.

First, do a fft of your signal so you know what the ‘valid’ frequency components are, and what constitutes noise. This is necessary for you to design your filters. I would begin with a Butterworth bandpass filter to filter out the d-c offset and baseline drift at the low end, and high-frequency noise at the high end. That should leave you with a relatively smooth signal. (My filter design procedure uses the Signal Processing Toolbox, but I include it here anyway for your reference: How to design a lowpass filter for ocean wave data in Matlab?.)

Then use the File Exchange PeakFinder on the negative of your signal to find the troughs, since that seems to be where you want to divide it. Use those indices on your original filtered signal (not the negative of it) to create a cell array of the segments of your signal (since it is unlikely all will be the same lengths).

If the baseline, baseline drift, trend or other aspects of your signal are important, you can specifically filter and retrieve those and save them separately.

[Image Analyst]

Star, what are the advantages of using Butterworth over sgolayfilt() (a local sliding polynomial filter)? Or just erasing high frequencies in the Fourier domain? What are the appropriate/best conditions for using each of these?

What are the advantages, if any, of using the File Exchange PeakFinder over the findpeaks() in the Signal Processing Toolbox?

[Star Strider]

‘but unfortunately a lot of them require the use of Toolboxes I do not have (Signal Processing, Symbolic, Image Processing...)... I use Matlab 2015b with Simulink and the Control System Toolbox.’

Short answer: He doesn’t have the Signal Processing Toolbox!