Remove trend and detect peaks in a photoplethysmogram(PPG) signal

There seems to be a sinusoidal like trend in my PPG signal. I believe this is due to respiration. I would like to detect the PPG peaks(maxima) in this signal after removing the trend or even better detect the peaks without removing the trends. Please help to provide a sample code as I am kind of new to signal processing for bio-signals.

clear all;
clc;

readData = csvread('p3_sit.csv',0,0);

ppg_head = readData(:,1).';

fs = 960; 
rec_time_minutes = (length(ppg_head)-1)/60000 %recording time calculation

wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,ppg_head);

N1 = length(ppg_head_data);
t1 = [0:N1-1]/fs;
%{comment here %}
figure(1)
plot(t1,ppg_head_data);
xlabel('second');ylabel('Volts');title('PPG Head Signal')

Answer 1

[Star Strider]

If you want to eliminate the baseline wander and offset, change to a highpass filter with a different cutoff frequency:

wn=5/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn,'high'); 
ppg_head_data = filter(b,a,ppg_head);

There are much better ways to design and implement filters. Here is one:

t = linspace(0, numel(ppg_head), numel(ppg_head))/fs;

Fs = 960;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = 24.5/Fn;                                               % Stopband Frequency (Normalised)
Ws = 25.0/Fn;                                               % Passband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                             % Filter Order
[z,p,k] = cheby2(n,Rs,Ws,'high');                           % Filter Design, Sepcify Bandstop
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability

figure(2)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot

ppg_head_data = filtfilt(sos, g, ppg_head);                 % Filter Signal

Use the Signal Processing Toolbox findpeaks (link) function on your filtered data to locate the peaks:

[pks,locs] = findpeaks(ppg_head_data, 'MinPeakHeight',2E+5);

You can use the ‘locs’ index vector to refer to the peaks and times in the original unfiltered signal.

[Kahuthaman Silvadorai]

Hi Star Strider. Thank you so much for the assistance. But in your code you have actually filtered out the desired signal instead and plotted the noise. This may be because I have mislead you from my previous code example. But my desired signal will be the PPG signal(that looks like rectified sinusoidal signal as in my attached plot) instead.

[Star Strider]

My pleasure.

It is still not clear what your filtered signal should look like.

You can always adjust the passband and stopband frequencies in my filter. In a highpass design, ‘Ws’ must be lower than ‘Wp’, and will pass all frequencies higher than ‘Ws’.

Also try changing 'high' to 'low' to create a lowpass filter from it if that does what you want. Note that in a lowpass design, ‘Ws’ must be higher than ‘Wp’, and the result will keep the constant offset and frequencies lower than ‘Ws’.

[Kahuthaman Silvadorai]

If you zoom in to my plot(in the initial code), you will notice some dark curvy peaks. I will need to filter them from the noise. I am unable to desgin the filter as the FFT of the signal is at 0hz. Not sure where to start. I did try ajusting your passband and stopband, but it doesnt seem to work for me. Please see the attachment to this answer.

[Star Strider]

This works:

Fs = 960;                                                   % Sampling Frequency (Hz)
t = linspace(0, numel(ppg_head), numel(ppg_head))/Fs;

Ts = mean(diff(t));                                         % Sampling Interval
Fs = 1/Ts;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
L = numel(t);                                               % Signal Length
ppg_head = ppg_head - mean(ppg_head);
FTsig = fft(ppg_head)/L;                                    % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector

figure(1)
plot(Fv, abs(FTsig(Iv))*2)
grid
% axis([0  3    ylim])

Fs = 960;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [1.5 10.00]/Fn;                                        % Passband Frequency (Normalised)
Ws = [1.0 10.50]/Fn;                                        % Stopband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                             % Filter Order
[z,p,k] = cheby2(n,Rs,Ws);                                  % Filter Design, Sepcify Bandpass
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability

figure(2)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot

ppg_head_data = filtfilt(sos, g, ppg_head);                 % Filter Signal

[pks,locs] = findpeaks(ppg_head_data, 'MinPeakHeight',1E+5, 'MinPeakDist',350);

figure(3)
plot(t, ppg_head_data, '-b')
hold on
plot(t(locs), pks, '^r')
hold off
axis([0  20.6    ylim])

You may need to experiment with the findpeaks arguments. I only looked at a section of the signal in detail.

[Kahuthaman Silvadorai]

Awesome. Thank you so much Star Strider. One more question for my learning sake, how did you compute the following values below. Can you elaborate on this: Wp = [1.5 10.00]/Fn; Ws = [1.0 10.50]/Fn; Rp = 1; Rs = 50;

[Star Strider]

As always, my pleasure.

I designed the filter first from the Fourier transform results, then refined it by simple experimentation. This is typical for biomedical signal processing.

[Kahuthaman Silvadorai]

oh ok..Just curios how to did come out with the stopband(1.0 10.50) values. Im not so good in interpreting FFT plots..And also do you have any recommendation for my other question: https://ch.mathworks.com/matlabcentral/answers/381643-unable-to-filter-noise-and-powerline-interference-from-head-and-chest-ecg-before-averaging

[Star Strider]

For a bandpass filter, the stopband frequencies must be ‘outside’ the passband frequencies (less than and greater than respectively). I usually choose a relatively small transition region, here 0.5 Hz, especially at low frequencies, to be certain that the stopbands do not include 0 Hz and produce a stable filter, also my reason for preferring Chebyshev Type II designs.

I read your other Question when you first posted it. I believe that is primarily an instrumentation problem. You need to have a ‘ground’ or reference electrode, usually the right leg lead, to eliminate powerline and other external signal interference. The internal differential amplifiers will then eliminate the common-source noise. This applies to EEG and EMG as well as EKG signal acquisition. All good biomedical instrumentation textbooks, such as John G. Webster’s, describe the theory and implementation of these techniques.

[Kahuthaman Silvadorai]

Hi..I did have a reference electrode placed at glabella. But the 50-60Hz noise seems to be unavoidable. Is it a way for me to eliminate these noise by filtering ? The recording instrument that I am using is Bioradio150.

[Star Strider]

One thing that you absolutely must do is to use shielded cables on all your electrode leads. These are small coaxial cables that allow you to ground the shield. I looked up BioRadio website, specifically the ECG system. I cannot find any specific information on the lead systems, other than it offers shielded cables for the ECG system. The records on the site do not display powerline noise.

It would most appropriate to address that problem with GLINT technical support.

I would do that before using filters, since a filter to ‘notch out’ the powerline frequency also notches out your signals at that frequency, and most physiological signals have some information at the powerline frequency.

[Syed Muhammad Abubakar]

Hi!

In this example you used IIR filter to remove the noise. I want to ask which FIR filter will depicts same results as of the Chebyshev Type II 2nd order filter used here? I want to reproduce the results of the first code shared by Star Strider.

Thanks in advance.

[632541]

Hi Star strider,

Why do you do subtract the mean signal from the signal itself?

ppg_head = ppg_head - mean(ppg_head);

What is the significance of it?

Can it be skipped?

Please answer these.

Thank you..

[Star Strider]

Syed Muhammad Abubakar — I woulld use kaiserord with fir1 to design a FIR filter. Use filtfilt to do the actual filtering.

Kaveri A — Subtracting the mean before calculating and ploting the fft makes the other peaks more easily visible. It is not necessary, however it makes the fft result easier to interpret.

.

[632541]

Anybody please answer this

[Star Strider]

Kaveri A — Why the desperation? I already did!

If that doesn’t resolve the issue, post a new Question!

[632541]

Thanks Star Strider...

But the frequency of the signal in fft plot chenges in my case without subtracting mean..

Will you please see the difference and tell me..

Subtracting

Without subtracting the mean

Thank you..

[632541]

No Star Strider ,actually I didn't refresh the page that's why didn't get update

[Star Strider]

Kaveri A — The frequency does not change, although it is difficult to see the blue line in the plot calculated without subtracting the mean. Note that the mean value (d-c offset) is about , while the next largest peak is aboout so the other peaks are essentially invisible.

Subtracting the mean before calculating the fft is easiest. Plotting the fft power in dB, or amplitude with semilogy, would show the peaks as well.

[632541]

Thanks a lot !

[632541]

Hi Star Strider ,

In the code ,

FTsig = fft(ppg_head)/L;

Can you please tell why divided by L?

What does dividing number of samples signify ?

And how is it related to signal amplitude?

[Star Strider]

@Kaveri A — It normalises it by the signal length. (It’s also discussed in the fft documentation.)

Answer 2

[Peter H Charlton]

Since the previous answer, a new toolbox has become available to detect peaks in a photoplethysomogram signal:

https://uk.mathworks.com/matlabcentral/fileexchange/107644-ppg-beats

After downloading the toolbox and adding it to the Matlab path (as described here), the toolbox can be used to detect peaks in this signal as follows (this is a slight modification of one of the tutorials in the toolbox documentation, here):

% setup using the code provided in the question
readData = csvread('p3_sit.csv',0,0);
ppg_head = readData(:,1).';
fs = 960; 
wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,ppg_head);

% format data as required by the toolbox
S.v = ppg_head_data(:);   % PPG samples in a column vector 
S.fs = fs;                % sampling frequency in Hz

% detect PPG beats
beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats

% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

This produces the following plot, which shows beats detected in the signal:

Zooming in to a 10-second segment:

xlim([10, 20]); h = findobj(gca,'Type','line'); set(h, 'LineWidth', 2)

This plot shows the detected peaks (red circles).

Disclaimer: I'm one of the authors of this PPG beat detection toolbox.

[Patri Valcárcel]

Good evening,

I'm trying to use your explanation for detecting PPG peaks in signal and I get an error 'Unrecognized function or variable 'detect_ppg_beats' ' where does that function comes from?

Thank you in advance.

[Peter H Charlton]

Thanks for your question.

Please note that the toolbox must be downloaded, and then added to the Matlab path, before using the above code.

I expect this error is due to the toolbox not being visible on the Matlab path (either because it hasn't been downloaded, or because the downloaded files haven't been added to the Matlab path).

I've now added instructions on how to download the files and add them to the Matlab path here. I'll also add these steps to my original answer.

[Patricia]

Thank you very much for the answer!

Excuse me, it is one of my first times using MATLAB and I am a bit lost. May I ask you something else? I have two signals, one coming from a camera (it is remote PPG, 32FPS) and another one coming from a conventional sensor (so it is typical PPG signal, fs=60) and I have to process them and find the peaks and IBIs on each signal. As I told you it is one of my first times working with these signals and I studied the toolbox to use it but I was wondering if the previous code would also work for my case. I have resampled (downsample) the PPG signal to the remote PPG signal frequency (60 to 32) to have them aligned and be able to compare them but once I try to use the lowpass filter you propose on the previous answers I see no diference on the output signal, moreover, I guess I need more pre-processing than what I am already doing since I think the algorithm is not detecting the peaks correctly (at least on the conventional ppg signal) because of my pre-processing. (And last question, can I use this code for both of the signals? I think it would work with both but I still have that question in mind)

Excuse me for the long question but I am willing to learn more and be able to finish this project! :)

I am attaching the code so far here and the data on the question, thank you in advance!

code:

%remote signal
Signals = load('RPPG.mat')
RPPGSignal = Signals.ProcessedData.RPPG;
fs1= 32;

%sensor signal
MonHearthRate=readtable('PPG_sensor.csv');
heartRate= table2array(MonHearthRate);
PPGSignal= heartRate.';
fs2=60;
%% resampling (downsample the signal coming from the sensor to rPPG (60Hz to 32HZ)
PPGSignal_resampled = (resample(PPGSignal,8,15)).'; %8/15 * 60Hz = 32Hz
PPGSignal_resampled = PPGSignal_resampled.';
%% Preparing the PPG signal
ppg_noOff = detrend(PPGSignal_resampled);


%% Detect peaks on noisy PPG Signals - Only on sensor PPG signal

% setup 
fs = 32; 
wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,y);

% format data as required by the toolbox
S.v = y(:);   % PPG samples in a column vector 
S.fs = fs;                % sampling frequency in Hz
% detect PPG beats
beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats
% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

[Peter H Charlton]

This code seems to work. Change 'do_remote' to zero to process the PPG sensor signal.

%% Load data

%remote signal
Signals = load('RPPG.mat');
RPPGSignal = Signals.ProcessedData.RPPG;
fs1= 32;

%sensor signal
MonHearthRate=readtable('PPG_sensor.csv');
heartRate= table2array(MonHearthRate);
PPGSignal= heartRate.';
fs2 = 60;

%% format data as required by the toolbox

do_remote = 1;
if do_remote
    S.v = RPPGSignal;
    S.fs = fs1;                % sampling frequency in Hz
else
    S.v = PPGSignal;
    S.fs = fs2;
end

%% detect PPG beats

beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats

%% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

[Patricia]

Thank you again! :) So should I not pre-process the signal? And is there any option on the toolbox to calculate the IBIs?

[Peter H Charlton]

I don't know whether pre-processing the signals would improve beat detection. In our recent article on PPG beat detection we pre-processed the signals by bandpass filtering between 0.67 and 8 Hz, which could help.

Inter-beat intervals (IBIs) can be calculated as follows:

ibis = diff(t(peaks));

[Patricia]

Hi again! the previous code gives me an error on the detect_ppg_beats as 'Too many input arguments' By what can this be caused?

[Peter H Charlton]

Good point - now corrected

[Patricia]

Thank you! :) And would you recommed to use the preprocess_ppg_signal function as a method to prepare the signal? Or just the bandpass filtering between 0.67 and 8 Hz?

[Patricia]

Hi again and excuse me for reopening the question. I was trying the performance assessment of the detector using the 'capnobase' option and once I get it working I do not get any result as a return, stops at this point:

any idea on what can I be doing wrong?

Thank you again.