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

141 views (last 30 days)
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')

Accepted Answer

Star Strider
Star Strider on 5 Feb 2018
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.
  21 Comments
Star Strider
Star Strider on 21 Apr 2021
@Kaveri A — It normalises it by the signal length. (It’s also discussed in the fft documentation.)

Sign in to comment.

More Answers (1)

Peter H Charlton
Peter H Charlton on 25 Nov 2022
Edited: Peter H Charlton on 28 Dec 2022 at 6:52
Since the previous answer, a new toolbox has become available to detect peaks in a photoplethysomogram signal:
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:
PPG signal with peaks detected
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.
  10 Comments
Patricia
Patricia on 20 Jan 2023 at 13:59
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!