Remove trend and detect peaks in a photoplethysmogram(PPG) signal
66 views (last 30 days)
Show older comments
Kahuthaman Silvadorai
on 5 Feb 2018
Commented: Patricia
on 20 Jan 2023
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')
0 Comments
Accepted Answer
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
632541
on 21 Apr 2021
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
on 21 Apr 2021
More Answers (1)
Peter H Charlton
on 25 Nov 2022
Edited: Peter H Charlton
on 28 Dec 2022
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:
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
on 12 Jan 2023
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
on 20 Jan 2023
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.
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!