SciCodePile/SciCode-Domain-Code
DATA1: Domain-Specific Code Dataset Dataset Overview DATA1 is a large-scale domain-specific code dataset focusing on code samples from interdisciplinary fields such as biology, chemistry, materials science, and related areas. The dataset is collected and organized from GitHub repositories, covering 178 different domain topics with over 1.1 billion lines of code. Dataset Statistics Total Datasets: 178 CSV files Total Data Size: ~115 GB Total Lines… See the full description on the dataset page: https://huggingface.co/datasets/SciCodePile/SciCode-Domain-Code.
42.4k
1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Biosensors","Cassey2016/PPG_Peak_Detection","main.m",".m","3040","58","% =========================================================================3% Below functions are the implementation for the comparison methods in4% paper:5% Han, Dong, Syed K. Bashar, Jesús Lázaro, Fahimeh Mohagheghian, 6% Andrew Peitzsch, Nishat Nishita, Eric Ding, Emily L. Dickson, 7% Danielle DiMezza, Jessica Scott, Cody Whitcomb, Timothy P. Fitzgibbons, 8% David D. McManus, and Ki H. Chon. 2022. 9% ""A Real-Time PPG Peak Detection Method for Accurate Determination of 10% Heart Rate during Sinus Rhythm and Cardiac Arrhythmia"" 11% Biosensors 12, no. 2: 82. https://doi.org/10.3390/bios12020082 12%13% Please cite our paper if you used our implementation code. Thank you.14% Author: Dong Han (dong.han@uconn.edu), 01/31/2022.15% =========================================================================16 17% -------------------------------------------------------------------------18% Input: 19% PPG_raw_buffer: should be 30-sec segment.20% fs_PPG_raw: the sampling frequency of the PPG_raw_buffer.21% -------------------------------------------------------------------------22%% Preparation of PPG signal:23addpath('.\func')24[PPG_buffer,fs_PPG] = my_func_prep_PPG_buffer(PPG_raw_buffer,fs_PPG_raw);25 26%% Method 1: implemented method 1-a27V_max_flag = true; % true == upper peak detection.28addpath('.\method_01_and_02');29output_upper_Shin_2009 = my_peak_compare_Shin_2009(PPG_buffer,fs_PPG,V_max_flag); % Implementation of Shin 2009 paper.30 31%% Method 2: implemented method 1-b32V_max_flag = false; % false == lower peak detection.33output_lower_Shin_2009 = my_peak_compare_Shin_2009(PPG_buffer,fs_PPG,V_max_flag); % Implementation of Shin 2009 paper.34 35%% Method 3 & 4: implemented method 2, it has two output peaks in ""output_Elgendi_1_2013""36delta = 0.5; % it was 0.1 as mentioned in the paper. But I think 0.5 works better (0.5 is in the billauer's website).37addpath('.\method_03_and_04');38[output_Elgendi_1_2013] = my_Elgendi_2013_method_I_peakdet(PPG_buffer, delta, fs_PPG);39 40%% Method 5: first derivative and adaptive thresholding method in Li et al. [4] and Elgendi's paper [3]41abpsig = resample(PPG_buffer,fs_abpsig,fs_PPG_buffer); % upsampling it to 125 Hz.42addpath('.\method_05');43[output_Elgendi_2_2013] = my_func_ppg_peakdet_method_05_Elgendi_2013_method_II(abpsig,fs_abpsig);44 45%% Method 6: implemented method 446fs_abp = 250; % Hz.47abp = resample(PPG_buffer,fs_abp,fs_PPG); % upsampling it to 125 Hz.48addpath('.\method_06');49[output_Elgendi_3_2013] = my_Elgendi_2013_method_III_peakdet(abp,fs_abp);50 51%% Method 7: event-related moving averages with dynamic threshold method in Elgendi et al.'s paper [3] 52addpath('.\method_07');53[output_Elgendi_4_2013] = my_func_ppg_peakdet_method_07_Elgendi_2013_method_IV(-PPG_raw_buffer,fs_PPG_raw);54 55%% Method 8 & 9: peak detection on Stationary Wavelet Transform of PPG signal56fs_swt = 125; % Hz.57PPG_swt = resample(PPG_buffer,fs_swt,fs_PPG); % upsampling it to 125 Hz.58addpath('.\method_08_and_09');59[output_Vadrevu_1_2019,output_Vadrevu_2_2019] = my_Vadrevu_2019_peakdet(PPG_swt,fs_swt);","MATLAB"
60"Biosensors","Cassey2016/PPG_Peak_Detection","method_07/my_func_ppg_peakdet_method_07_Elgendi_2013_method_IV.m",".m","6505","156","function output_Elgendi_4_2013 = my_func_ppg_peakdet_method_07_Elgendi_2013_method_IV(raw_PPG,fs_PPG)61% =========================================================================62% This is my implementation of the method IV in this paper:63% Elgendi, Mohamed, et al. 64% ""Systolic peak detection in acceleration photoplethysmograms measured from 65% emergency responders in tropical conditions."" PLoS One 8.10 (2013): e76585.66%67% Implemented by Dong Han on 03/02/2020.68%69% Please cite our paper if you used this code:70% Han, Dong, Syed K. Bashar, Jesús Lázaro, Fahimeh Mohagheghian, 71% Andrew Peitzsch, Nishat Nishita, Eric Ding, Emily L. Dickson, 72% Danielle DiMezza, Jessica Scott, Cody Whitcomb, Timothy P. Fitzgibbons, 73% David D. McManus, and Ki H. Chon. 2022. 74% ""A Real-Time PPG Peak Detection Method for Accurate Determination of 75% Heart Rate during Sinus Rhythm and Cardiac Arrhythmia"" 76% Biosensors 12, no. 2: 82. https://doi.org/10.3390/bios12020082 77%78% Please cite our paper if you used our code. Thank you.79% =========================================================================80%% pre-processing - bandpass filtering81 [b, a] = butter(2,[0.5 8]/(fs_PPG/2)); % 2nd order bandpass filter 0.5-8Hz;82 filtered_PPG = filtfilt(b, a, raw_PPG); % zero-phase filter.83 filtered_PPG = filtered_PPG ./ std(filtered_PPG); % normalizing data is very important for my peak detection.84 filtered_PPG = filtered_PPG - mean(filtered_PPG);85 86 debugging_plot_flag = false; % only for plotting debugging figures. 87 88 % clip the signal by keeping the signal above zero.89 % I do not want to do this, so i will move all signal above zero.90 S_n = filtered_PPG;91% ---- Not following the paper to clip signal but move all signal above zero:92% if min(S_n) < 093% Z_n = S_n - min(S_n); % elevate signal above zero.94% else95% % the minimum of S_n is still above zero, so do nothing.96% Z_n = S_n;97% end98% ---- Following the paper: only keep the positive value:99 Z_n = S_n;100 Z_n(Z_n < 0) = 0; 101%% pre-processing - squaring102 y_n = (Z_n).^2; % element-wise power.103%% feature extraction - generating potential blocks using two moving averages104 W_1 = round(0.111 * fs_PPG); % mentioned as the paper by brute-force search.105 % first moving average:106% MA_peak = y_n; % for the beginning and ending signal, use the original signal.107% for nn = 1+round(W_1/2):length(raw_PPG)-round(W_1/2)108% temp_range = (nn-round(W_1/2)):(nn+round(W_1/2));109% MA_peak(nn) = sum(y_n(temp_range))/W_1;110% end111 MA_peak = movmean(y_n,W_1);112 113 % second moving average: 114 W_2 = round(0.667 * fs_PPG);115% MA_beat = y_n;116% for nn = 1+round(W_2/2):length(raw_PPG)-round(W_2/2)117% temp_range = (nn-round(W_2/2)):(nn+round(W_2/2));118% MA_beat(nn) = sum(y_n(temp_range))/W_2;119% end120 MA_beat = movmean(y_n,W_2);121%% classification - thresholding122 beta = 0.02; % from the paper, by brute force search.123 z_bar = mean(y_n);124 alpha = beta * z_bar; % offset level.125 THR_1 = MA_beat + alpha;126 127 Blocks_Of_Interest = zeros(size(MA_peak)); % I initial it as zero.128 for nn = 1:length(MA_peak)129 if MA_peak(nn) > THR_1(nn) % I think it is THR_1(nn).130 Blocks_Of_Interest(nn) = 0.1;131 else132 % since I inital block of interest as zero, so I do not need to133 % assign zero again.134 end135 end136 137 % searh for onset and offset of each block.138 count_blocks = 0;139 block_onset = NaN(size(MA_peak));140 block_offset = NaN(size(MA_peak));141 if any(Blocks_Of_Interest > 0) % there is a block exist.142 for nn = 1:length(MA_peak)143 if nn == 1 && Blocks_Of_Interest(nn) > 0144 % the first point is a block;145 count_blocks = count_blocks + 1; % since the block start from zero, I have to add the counter first.146 block_onset(count_blocks,1) = nn;147 elseif nn == length(MA_peak) && Blocks_Of_Interest(nn) > 0148 % end with a block:149 % no need to add count_blocks;150 block_offset(count_blocks,1) = nn;151 else152 if nn > 1153 if Blocks_Of_Interest(nn-1) == 0 && Blocks_Of_Interest(nn) > 0 % a jump means a new block.154 count_blocks = count_blocks + 1;155 block_onset(count_blocks,1) = nn;156 elseif Blocks_Of_Interest(nn-1) > 0 && Blocks_Of_Interest(nn) == 0 % a drop means the end of previous block.157 block_offset(count_blocks,1) = nn;158 end159 end160 end161 end162 else163 % there is no block existed. Check why.164% keyboard;165 HR_Elgendi_4_2013 = 0; % there is no peak location.166 S_peaks = 1;167 output_Elgendi_4_2013 = struct('filtered_PPG_Elgendi_4_2013',S_n,...168 'PPG_peak_loc_Elgendi_4_2013',S_peaks,...169 'HR_Elgendi_4_2013',HR_Elgendi_4_2013);170 return171 end172 173 block_onset(isnan(block_onset)) = []; % remove extra elements.174 block_offset(isnan(block_offset)) = []; % remove extra elements.175 if size(block_onset,1) ~= size(block_offset,1)176 % not same number of onset and offset, check here.177 keyboard;178 end179 180 if size(block_onset,1) ~= count_blocks181 keyboard;182 end183 S_peaks = NaN(count_blocks,1);184 THR_2 = W_1;185 186 for jj = 1:count_blocks187 block_idx = [block_onset(jj,1):block_offset(jj,1)];188 [~,I] = max(y_n(block_idx));189 S_peaks(jj,1) = block_onset(jj,1) + I - 1;190 end191 192 if debugging_plot_flag193 figure;194 plot(filtered_PPG);hold on;195 plot(S_peaks,y_n(S_peaks),'r.','markersize',10);196 plot(y_n);197 plot(MA_peak,'k:');198 plot(MA_beat,'r--');199 plot(THR_1,'g.-');200 plot(Blocks_Of_Interest*max(y_n)*10,'color',[0.5,0.5,0.5]); % grey color. I want to make block more obvious.201 202 legend('filtered PPG','peaks', 'squared PPG with clip to zero', 'MA peak', 'MA beat','THR 1', 'Blocks of Interest');203 end204 205 if isempty(S_peaks)206 HR_Elgendi_4_2013 = 0; % there is no peak location.207 S_peaks = 1;208 else209 HR_Elgendi_4_2013 = 60 * fs_PPG ./ diff(S_peaks); % calculate the HR.210 end211 212 output_Elgendi_4_2013 = struct('filtered_PPG_Elgendi_4_2013',S_n,...213 'PPG_peak_loc_Elgendi_4_2013',S_peaks,...214 'HR_Elgendi_4_2013',HR_Elgendi_4_2013);215end","MATLAB"
216"Biosensors","Cassey2016/PPG_Peak_Detection","method_01_and_02/my_peak_compare_Shin_2009.m",".m","21523","388","function [output_Shin_2009] = my_peak_compare_Shin_2009(raw_PPG,fs_PPG,V_max_flag)217% =========================================================================218% This function is the implementation of this paper:219% Shin, Hang Sik, Chungkeun Lee, and Myoungho Lee. 220% ""Adaptive threshold method for the peak detection of 221% photoplethysmographic waveform."" 222% Computers in biology and medicine 223% 39.12 (2009): 1145-1152.224% 225% Implemented by: Dong Han, on 02/10/2020.226%227% Please cite our paper if you used this code:228% Han, Dong, Syed K. Bashar, Jesús Lázaro, Fahimeh Mohagheghian, 229% Andrew Peitzsch, Nishat Nishita, Eric Ding, Emily L. Dickson, 230% Danielle DiMezza, Jessica Scott, Cody Whitcomb, Timothy P. Fitzgibbons, 231% David D. McManus, and Ki H. Chon. 2022. 232% ""A Real-Time PPG Peak Detection Method for Accurate Determination of 233% Heart Rate during Sinus Rhythm and Cardiac Arrhythmia"" 234% Biosensors 12, no. 2: 82. https://doi.org/10.3390/bios12020082 235%236% Please cite our paper if you used our code. Thank you.237% =========================================================================238 debugging_plot_flag = false; % debugging plot. Can be false if don't want to plot anything.239%% Section 2.4 PPG frequency analysis and filtering.240 241 % (1): high pass >= 0.5 Hz.242 [b, a] = butter(6,[0.5 20]/(fs_PPG/2)); % bandpass filter 0.5-10Hz, changed from 0.5-20 to 0.5-9 Hz at 11/21/2018243 raw_PPG = filtfilt(b, a, raw_PPG); % -> AC component244 raw_PPG = raw_PPG ./ std(raw_PPG); % normalizing data is very important for my peak detection.245 raw_PPG = raw_PPG - mean(raw_PPG);246%% Section 2.5 & 2.6 Peak detection algorithm & Adaptive threshold detection247 248 % (1): bandpass filtering, no moving average filter or wavelet249 % decomposition.250 filtered_PPG = raw_PPG;251 Fs = fs_PPG;252 253 % % ===== interpolation to 1kHz of PPG: =====254 % x = 1:length(filtered_PPG);255 % v = filtered_PPG;256 % 257 % upsample_Fs = 250;258 % xq = 1:Fs/upsample_Fs:length(filtered_PPG);259 % vq1 = interp1(x,v,xq);260 % 261 % filtered_PPG = vq1;262 % Fs = upsample_Fs; % upsampled to 1000 Hz.263 264 % figure265 % plot(x,v,'o',xq,vq1,':.');266 % xlim([0 max(xq)]);267 % title('(Default) Linear Interpolation');268 269 % (2): V_max270 % slope_k: k-th slope amplitude;271 % s_r: slope changing rate (empirically: V_max = -0.6);272 % V_n_1: previous peak amplitude;273 % std_PPG: standard deviation of entire PPG signal;274 % Fs: sampling frequency.275 276 filtered_PPG = filtered_PPG(:);277 slope_k = NaN(size(filtered_PPG)); % should be a column vector.278 peak_loc = NaN(size(filtered_PPG)); % the array to store PPG peak index.279 pk_idx = 1; % the counter of peaks.280 %% Section 2.7: Peak Correction281 refractory_period = 0.6 * Fs; % sec * sampling frequency, initial refractory period is 0.6 sec.282 283 temp_win_left = round(0.15 * Fs); % sec * sampling frequency. This is the search region for local minima or maxima detection. chose 0.15 sec because 0.3 sec == 200 BPM.284 temp_win_right = round(0.15 * Fs);285 286 287 if V_max_flag % doing upper peak detection.288 s_r = -0.6; 289 else290 s_r = 0.6;%0.6; % not positive because my signal is zero mean.291 % I need to make all bottom signal positive, so I am moving them up.292 % move_filter_amp = min(filtered_PPG) * (-1);293 % filtered_PPG = filtered_PPG + move_filter_amp + std(raw_PPG); % move the lowest value more than zero.294 end295 296 297 slope_meet_PPG_flag = false; % mark if the slope meet PPG.298 slope_lower_PPG_flag = false; % mark if slope is lower than PPG, once PPG amp is lower than slope, mark it back.299 prev_slope = NaN; % First, I want to test not decreasing with PPG amplitude.300 if debugging_plot_flag % debugging plot301 figure;302 plot(filtered_PPG);303 hold on;304 end305 for kk = 1:length(filtered_PPG)306 % this is for debugging:307 if kk == 2308 my_stop = 1;309 end310 if kk == 1 % initial the slope value311 if V_max_flag312 slope_k(1,1) = 0.2 * max(filtered_PPG);313 std_PPG = std(filtered_PPG);314 else315 slope_k(1,1) = 0.2 * min(filtered_PPG); % since my signal is zero mean, I start from the negative amp. % I added what I moved.316 std_PPG = -std(filtered_PPG);317 end318 % std_PPG = std(filtered_PPG);319 V_n_1 = slope_k(1,1);320 else321 if slope_meet_PPG_flag % slope has met PPG before.322 slope_k(kk,1) = filtered_PPG(kk,1);323 if V_max_flag % upper peak detection.324 if kk < 2 % in the second point of signal325 turn_point_flag = (slope_k(kk,1) < slope_k(kk-1,1)); % we met local maximum.326 else327 turn_point_flag = (slope_k(kk,1) < slope_k(kk-1,1)) & (slope_k(kk - 1,1) > slope_k(kk-2,1)); % we met local maximum.328 end329 else330 if kk < 2 % in the second point of signal331 turn_point_flag = (slope_k(kk,1) > slope_k(kk-1,1)); % we met local minimum.332 else333 turn_point_flag = (slope_k(kk,1) > slope_k(kk-1,1)) & (slope_k(kk - 1,1) < slope_k(kk-2,1)); % we met local minimum.334 end335 end336 337 if turn_point_flag % there is a turning point.338 if pk_idx > 1 % not the first peak339 % check local maxima or minima:340 if (kk - temp_win_left) < 1341 temp_left = 1;342 else343 temp_left = kk - temp_win_left;344 end345 346 if (kk + temp_win_right) > length(filtered_PPG)347 temp_right = length(filtered_PPG);348 else349 temp_right = kk + temp_win_right;350 end351 temp_win = temp_left:temp_right;352 local_m_check = filtered_PPG(temp_win);353 if V_max_flag354 temp_m_idx = find(local_m_check > slope_k(kk - 1,1)); % check if there is another maximum than detected, remember use k-1.355 else356 temp_m_idx = find(local_m_check < slope_k(kk - 1,1)); % check if there is another minimum than detected357 end358 359 if isempty(temp_m_idx) % there is no more max or min than this peak360 if (kk - peak_loc(pk_idx-1,1) > refractory_period) % it is not the first peak, and the second peak is outside refractory period. It should be kk, because I have not assign the peak to the array.361 peak_loc(pk_idx,1) = kk-1;362 V_n_1 = filtered_PPG(peak_loc(pk_idx-1,1),1);% previous peak amplitude %slope_k(kk-1,1);363 % update refractory period:364 refractory_period = 0.6 * (kk - peak_loc(pk_idx-1,1)); % current index minus peak location. update the refractory peroid before updating the peak counting.365 pk_idx = pk_idx + 1;366 367 % reset slope meet flag:368 slope_meet_PPG_flag = false;369 slope_k(kk,1) = slope_k(kk - 1,1) + s_r * ((V_n_1 + std_PPG) / Fs); 370 371 % ---- for checking lower slope -------372 temp_slope_check = s_r * ((V_n_1 + std_PPG) / Fs);373 if V_max_flag374 if temp_slope_check > 0 % upper peaks should be decreasing with negative slope.375 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;376 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;377 end378 else379 if temp_slope_check < 0 % upper peaks should be decreasing with negative slope.380 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;381 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;382 end383 end384 % -------------------------------------------385 if V_max_flag386 temp_slope_below_PPG_flag = slope_k(kk,1) < filtered_PPG(kk,1); % upper peak detection, so slope below signal.387 else388 temp_slope_below_PPG_flag = slope_k(kk,1) > filtered_PPG(kk,1); % lower peak detection, so slope above signal.389 end390 if temp_slope_below_PPG_flag % if slope is below PPG signal, we will reset slope value to PPG amplitude.391 slope_lower_PPG_flag = true; % slope is lower than PPG signal.392 prev_slope = slope_k(kk,1); % store the slope value now.393 slope_k(kk,1) = filtered_PPG(kk,1);394 end395 if debugging_plot_flag % debugging plot396 plot(kk,slope_k(kk,1),'r.');397 end398 399 else400 if (kk - peak_loc(pk_idx-1,1) <= refractory_period) % it is because of the refractory period that cause the no peak. It should be kk, because I have not assign the peak to the array. 401 slope_k(kk,1) = filtered_PPG(kk,1);% from the fig.3(c) in the paper, I see they are using the signal amplitude, not slope. 402 % no need to reset slope meet flag, waiting for403 % next turning point.404 if debugging_plot_flag % debugging plot405 plot(kk,slope_k(kk,1),'r.');406 end407 end408 end409 else % there are more peaks higher then current kk peak.410 if debugging_plot_flag % debugging plot411 plot(kk,slope_k(kk,1),'r.');412 end413 end414 else % the first peak, no need to check refractory period.415 % check local maxima or minima:416 if (kk - temp_win_left) < 1417 temp_left = 1;418 else419 temp_left = kk - temp_win_left;420 end421 422 if (kk + temp_win_right) > length(filtered_PPG)423 temp_right = length(filtered_PPG);424 else425 temp_right = kk + temp_win_right;426 end427 temp_win = temp_left:temp_right;428 local_m_check = filtered_PPG(temp_win);429 if V_max_flag430 temp_m_idx = find(local_m_check > slope_k(kk-1,1)); % check if there is another maximum than detected, always detect previous peak.431 else432 temp_m_idx = find(local_m_check < slope_k(kk-1,1)); % check if there is another minimum than detected433 end434 435 if isempty(temp_m_idx)436 peak_loc(pk_idx,1) = kk-1;437 if pk_idx > 1438 V_n_1 = filtered_PPG(peak_loc(pk_idx-1,1),1);439 else440 V_n_1 = slope_k(kk-1,1);% previous peak amplitude %slope_k(kk-1,1);441 end442 pk_idx = pk_idx + 1;443 444 % reset slope meet flag:445 slope_meet_PPG_flag = false;446 slope_k(kk,1) = slope_k(kk - 1,1) + s_r * ((V_n_1 + std_PPG) / Fs); 447 % ---- for checking lower slope -------448 temp_slope_check = s_r * ((V_n_1 + std_PPG) / Fs);449 if V_max_flag450 if temp_slope_check > 0 % upper peaks should be decreasing with negative slope.451 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;452 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;453 end454 else455 if temp_slope_check < 0 % upper peaks should be decreasing with negative slope.456 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;457 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;458 end459 end460 % -------------------------------------------461 if V_max_flag462 temp_slope_below_PPG_flag = slope_k(kk,1) < filtered_PPG(kk,1); % upper peak detection, so slope below signal.463 else464 temp_slope_below_PPG_flag = slope_k(kk,1) > filtered_PPG(kk,1); % lower peak detection, so slope above signal.465 end466 467 if temp_slope_below_PPG_flag % if slope is below PPG signal, we will reset slope value to PPG amplitude.468 slope_k(kk,1) = filtered_PPG(kk,1);469 end470 if debugging_plot_flag % debugging plot471 plot(kk,slope_k(kk,1),'r.');472 end473 else % there are more peaks higher then current kk peak.474 if debugging_plot_flag % debugging plot475 plot(kk,slope_k(kk,1),'r.');476 end477 end 478 % no need to calculate refractory period, because there is only one peak, at least two peaks can give this correctly:479 end480 else481 % turning point did not meet, so keep decreasing or482 % increasing the slope.483 % slope_k(kk,1) = slope_k(kk - 1,1) + s_r * ((V_n_1 + std_PPG) / Fs); 484 if debugging_plot_flag % debugging plot485 plot(kk,slope_k(kk,1),'r.');486 end487 end488 else % slope has not met PPG before. Keep decresing or increasing according to 'V_max_flag'.489 % if slope_lower_PPG_flag % if there is a slope lower than PPG before:490 % slope_k(kk,1) = prev_slope;491 % else492 slope_k(kk,1) = slope_k(kk - 1,1) + s_r * ((V_n_1 + std_PPG) / Fs);493 % ---- for checking lower slope -------494 temp_slope_check = s_r * ((V_n_1 + std_PPG) / Fs);495 if V_max_flag496 if temp_slope_check > 0 % upper peaks should be decreasing with negative slope.497 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;498 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;499 end500 else501 if temp_slope_check < 0 % upper peaks should be decreasing with negative slope.502 temp_slope_check = -s_r * ((V_n_1 + std_PPG) / Fs);%-temp_slope_check;503 slope_k(kk,1) = slope_k(kk - 1,1) + temp_slope_check;504 end505 end506 % -------------------------------------------507 508 % end509 % if slope_k(kk,1) < filtered_PPG(kk,1) % if slope is below PPG signal, we will reset slope value to PPG amplitude.510 % slope_lower_PPG_flag = true; % slope is lower than PPG signal.511 % prev_slope = slope_k(kk,1); % store the slope value now.512 % slope_k(kk,1) = filtered_PPG(kk,1);513 % elseif slope_k(kk,1) > filtered_PPG(kk,1) % slope is higher.514 % slope_lower_PPG_flag = false;515 % prev_slope = NaN; % reset the prev value.516 517 % end518 519 % if slope_lower_PPG_flag ~= 1 % if slope was not lower than PPG.520 % % -------------- Check if two lines will meet -----------------521 % PPG_x1 = kk - 1;522 % PPG_x2 = kk;523 % PPG_y1 = filtered_PPG(kk-1,1);524 % PPG_y2 = filtered_PPG(kk,1);525 % slope = s_r;526 % slope_y2 = slope_k(kk,1);527 % slope_y1 = slope_k(kk-1,1);528 % [meet_x] = my_slope_meet_PPG(PPG_x1,PPG_x2,PPG_y1,PPG_y2,slope,slope_y2,slope_y1);529 % 530 % slope_meet_PPG_flag = (ceil(meet_x) == kk);%(slope_k(kk,1) - filtered_PPG(kk,1)) < 0.1; % 0.3 is a testing value. %slope_k(kk,1) == filtered_PPG(kk,1) % slope meets the PPG signal.531 % end532 if V_max_flag533 slope_meet_PPG_flag = ((slope_k(kk,1) < filtered_PPG(kk,1)) & slope_k(kk - 1,1) > filtered_PPG(kk - 1,1));534 else535 slope_meet_PPG_flag = ((slope_k(kk,1) > filtered_PPG(kk,1)) & slope_k(kk - 1,1) < filtered_PPG(kk - 1,1)); % lower peak use inverse amplitude.536 end537 % -------------------------------------------------------------538 % I found I cannot use equal, because the PPG sampling539 % frequency is not so high.540 if slope_meet_PPG_flag541 slope_k(kk,1) = filtered_PPG(kk,1); % starts from the next index, slope == PPG amplitude.542 else543 % don't need to do anything.544 if slope_lower_PPG_flag ~= 1 % there was no slope lower than PPG before.545 if V_max_flag546 slope_lower_PPG_flag = ((slope_k(kk,1) < filtered_PPG(kk,1)) & slope_k(kk - 1,1) == filtered_PPG(kk - 1,1)); % beginning part has same amplitude, but the ending part slope is lower.547 else548 slope_lower_PPG_flag = ((slope_k(kk,1) > filtered_PPG(kk,1)) & slope_k(kk - 1,1) == filtered_PPG(kk - 1,1)); % lower peak use inverse amplitude.549 end550 if slope_lower_PPG_flag551 prev_slope = slope_k(kk,1); % store the slope value now.552 slope_k(kk,1) = filtered_PPG(kk,1); % starts from the next index, slope == PPG amplitude.553 end554 else % there was slope lower than PPG before.555 556 if V_max_flag557 temp_PPG_below_slope_flag = filtered_PPG(kk,1) < prev_slope; % upper peak detection, so PPG below slope.558 else559 temp_PPG_below_slope_flag = filtered_PPG(kk,1) > prev_slope; % lower peak detection, so PPG above slope.560 end561 562 if temp_PPG_below_slope_flag % PPG is lower than prev slope.563 slope_k(kk,1) = prev_slope; % stop tracking PPG amp.564 slope_lower_PPG_flag = false; % reset the lower PPG flag.565 prev_slope = NaN;566 else567 slope_k(kk,1) = filtered_PPG(kk,1); % keep tracking PPG amp.568 end569 end570 end571 if debugging_plot_flag % debugging plot572 plot(kk,slope_k(kk,1),'r.');573 end574 end575 end576 577 end578 % ================== IMPORTANT: clean up NaN value ========================579 peak_loc(isnan(peak_loc)) = []; % remove empty peak loc.580 if V_max_flag % doing upper peak detection.581 582 else583 % moving signal back.584 % filtered_PPG = filtered_PPG - move_filter_amp - std(raw_PPG); % move the lowest value more than zero.585 % slope_k = slope_k - move_filter_amp - std(raw_PPG); % move the slope as well.586 end587 588 if debugging_plot_flag % debugging plot589 plot(peak_loc,filtered_PPG(peak_loc),'ko');590 end591 592 if isempty(peak_loc)593 HR_Shin_2009 = 0; % there is no peak location.594 peak_loc = 1;595 else596 HR_Shin_2009 = 60 * Fs ./ diff(peak_loc); % calculate the HR.597 end598 599 output_Shin_2009 = struct('PPG_peak_loc_Shin_2009',peak_loc,...600 'slope_Shin_2009',slope_k,...601 'filtered_PPG_Shin_2009',filtered_PPG,...602 'HR_Shin_2009',HR_Shin_2009); 603end","MATLAB"
604"Biosensors","Cassey2016/PPG_Peak_Detection","method_05/my_func_ppg_peakdet_method_05_Elgendi_2013_method_II.m",".m","11993","412","function [output_Elgendi_2_2013] = my_func_ppg_peakdet_method_05_Elgendi_2013_method_II(raw_PPG,fs_PPG) 605% -------------------------------------------------------------------------606% This peak detection function was mentioned in this paper:607% Elgendi, Mohamed, et al. 608% ""Systolic peak detection in acceleration photoplethysmograms measured from 609% emergency responders in tropical conditions."" PLoS One 8.10 (2013): e76585.610% 611 [onsetp,peakp,dicron,abpsig] = delineator(raw_PPG,fs_PPG);612% -------------------------------------------------------------------------613 614 if isempty(peakp) % there is no peak detected:615 HR_Elgendi_2_2013 = 0; % there is no peak location.616 peakp = 1;617 else618 HR_Elgendi_2_2013 = 60 * fs_PPG ./ diff(peakp); % calculate the HR.619 end620 621 622 output_Elgendi_2_2013 = struct('PPG_peak_loc_Elgendi_2_2013',peakp,...623 'HR_Elgendi_2_2013',HR_Elgendi_2_2013,...624 'filtered_PPG_Elgendi_2_2013',abpsig);625end626 627function [onsetp,peakp,dicron,abpsig] = delineator(abpsig,abpfreq)628% Below was copied from Mathwords File Exchange ""Pulse Waveform Delineator"": 629% https://www.mathworks.com/matlabcentral/fileexchange/29484-pulse-waveform-delineator630 631% This program is intended to delineate the fiducial points of pulse waveforms632% Inputs:633% abpsig: input as original pulse wave signals;634% abpfreq: input as the sampling frequency;635% Outputs:636% onsetp: output fiducial points as the beginning of each beat;637% peakp: output fiducial points as systolic peaks;638% dicron: output fiducial points as dicrotic notches;639 640% Its delineation is based on the self-adaptation in pulse waveforms, but641% not in the differentials.642 643% Reference:644% BN Li, MC Dong & MI Vai (2010) 645% On an automatic delineator for arterial blood pressure waveforms646% Biomedical Signal Processing and Control 5(1) 76-81.647 648% LI Bing Nan @ University of Macau, Feb 2007649% Revision 2.0.5, Apr 2009650 651%Initialization652peakIndex=0;653onsetIndex=0;654dicroIndex=0;655stepWin=2*abpfreq;656closeWin=floor(0.1*abpfreq); %invalide for pulse beat > 200BPM657 658sigLen=length(abpsig);659 660peakp=[];661onsetp=[];662dicron=[];663 664%lowpass filter at first665coh=25; %cutoff frequency is 25Hz666coh=coh*2/abpfreq;667od=3; %3rd order bessel filter668[B,A]=besself(od,coh);669abpsig=filter(B,A,abpsig);670abpsig=10*abpsig;671 672abpsig=smooth(abpsig);673 674%Compute differentials675ttp=diff(abpsig);676diff1(2:sigLen)=ttp;677diff1(1)=diff1(2);678diff1=100*diff1;679clear ttp;680diff1=smooth(diff1);681 682if sigLen>12*abpfreq683 tk=10;684elseif sigLen>7*abpfreq685 tk=5;686elseif sigLen>4*abpfreq687 tk=2;688else689 tk=1;690end691 692%Seek avaerage threshold in original signal693if tk>1 %self-learning threshold with interval sampling694 tatom=floor(sigLen/(tk+2));695 for ji=1:tk %search the slopes of abp waveforms696 sigIndex=ji*tatom;697 tempIndex=sigIndex+abpfreq;698 [tempMin,jk,tempMax,jl]=seeklocales(abpsig,sigIndex,tempIndex);699 tempTH(ji)=tempMax-tempMin;700 end701 abpMaxTH=mean(tempTH);702else703 [tempMin,jk,tempMax,jl]=seeklocales(abpsig,closeWin,sigLen);704 abpMaxTH=tempMax-tempMin;705end706clear j*;707clear t*;708 709abpMaxLT=0.4*abpMaxTH;710 711%Seek pulse beats by MinMax method712% diffIndex=1;713diffIndex=closeWin; %Avoid filter distortion714 715while diffIndex<sigLen716 tempMin=abpsig(diffIndex); %Initialization717 tempMax=abpsig(diffIndex);718 tempIndex=diffIndex;719 tpeakp=diffIndex; %Avoid initial error720 tonsetp=diffIndex; %Avoid initial error721 722 while tempIndex<sigLen723 %If no pulses within 2s, then adjust threshold and retry724 if (tempIndex-diffIndex)>stepWin725% tempIndex=diffIndex-closeWin;726 tempIndex=diffIndex;727 abpMaxTH=0.6*abpMaxTH;728 if abpMaxTH<=abpMaxLT729 abpMaxTH=2.5*abpMaxLT;730 end731 break;732 end733 734 if (diff1(tempIndex-1)*diff1(tempIndex+1))<=0 %Candidate fiducial points735 if (tempIndex+5)<=sigLen736 jk=tempIndex+5;737 else738 jk=sigLen;739 end740 if (tempIndex-5)>=1741 jj=tempIndex-5;742 else743 jj=1;744 end745 746 %Artifacts of oversaturated or signal loss?747 if (jk-tempIndex)>=5748 for ttk=tempIndex:jk749 if diff1(ttk)~=0750 break;751 end752 end753 if ttk==jk754 break; %Confirm artifacts755 end756 end757 758 if diff1(jj)<0 %Candidate onset759 if diff1(jk)>0760 [tempMini,tmin,ta,tb]=seeklocales(abpsig,jj,jk);761 if abs(tmin-tempIndex)<=2762 tempMin=tempMini;763 tonsetp=tmin;764 end765 end766 elseif diff1(jj)>0 %Candidate peak767 if diff1(jk)<0768 [tc,td,tempMaxi,tmax]=seeklocales(abpsig,jj,jk);769 if abs(tmax-tempIndex)<=2770 tempMax=tempMaxi;771 tpeakp=tmax;772 end773 end774 end775 776 if ((tempMax-tempMin)>0.4*abpMaxTH) %evaluation777 if ((tempMax-tempMin)<2*abpMaxTH)778 if tpeakp>tonsetp779 %If more zero-crossing points, further refine!780 ttempMin=abpsig(tonsetp);781 ttonsetp=tonsetp;782 for ttk=tpeakp:-1:(tonsetp+1)783 if abpsig(ttk)<ttempMin784 ttempMin=abpsig(ttk);785 ttonsetp=ttk;786 end787 end788 tempMin=ttempMin;789 tonsetp=ttonsetp;790 791 if peakIndex>0792 %If pulse period less than eyeclose, then artifact793 if (tonsetp-peakp(peakIndex))<(3*closeWin)794 %too many fiducial points, then reset795 tempIndex=diffIndex; 796 abpMaxTH=2.5*abpMaxLT;797 break;798 end799 800 %If pulse period bigger than 2s, then artifact801 if (tpeakp-peakp(peakIndex))>stepWin802 peakIndex=peakIndex-1;803 onsetIndex=onsetIndex-1;804 if dicroIndex>0805 dicroIndex=dicroIndex-1;806 end807 end808 809 if peakIndex>0810 %new pulse beat811 peakIndex=peakIndex+1;812 peakp(peakIndex)=tpeakp;813 onsetIndex=onsetIndex+1;814 onsetp(onsetIndex)=tonsetp;815 816 tf=onsetp(peakIndex)-onsetp(peakIndex-1);817 818 to=floor(abpfreq./20); %50ms819 tff=floor(0.1*tf);820 if tff<to821 to=tff;822 end823 to=peakp(peakIndex-1)+to;824 825 te=floor(abpfreq./2); %500ms826 tff=floor(0.5*tf);827 if tff<te828 te=tff;829 end830 te=peakp(peakIndex-1)+te;831 % Dong added on 05/07/2020:832 % For MIMIC III PACPVC 3_2, ii = 25.833 if te > length(diff1)834 te = length(diff1);835 end836 tff=seekdicrotic(diff1(to:te));837 if tff==0838 tff=te-peakp(peakIndex-1);839 tff=floor(tff/3);840 end841 dicroIndex=dicroIndex+1;842 dicron(dicroIndex)=to+tff;843 844 tempIndex=tempIndex+closeWin;845 break;846 end847 end848 849 if peakIndex==0 %new pulse beat850 peakIndex=peakIndex+1;851 peakp(peakIndex)=tpeakp;852 onsetIndex=onsetIndex+1;853 onsetp(onsetIndex)=tonsetp;854 855 tempIndex=tempIndex+closeWin;856 break;857 end858 end859 end860 end861 end862 863 tempIndex=tempIndex+1; %step forward864 end865 866% diffIndex=tempIndex+closeWin; %for a new beat867 diffIndex=tempIndex+1;868end869 870if isempty(peakp),return;end871%Compensate the offsets of lowpass filter872sigLen=length(peakp);873for diffIndex=1:sigLen %avoid edge effect874 tempp(diffIndex)=peakp(diffIndex)-od;875end876ttk=tempp(1);877if ttk<=0878 tempp(1)=1;879end 880clear peakp;881peakp=tempp;882clear tempp;883 884sigLen=length(onsetp);885for diffIndex=1:sigLen886 tempp(diffIndex)=onsetp(diffIndex)-od;887end888ttk=tempp(1);889if ttk<=0890 tempp(1)=1;891end 892clear onsetp;893onsetp=tempp;894clear tempp;895 896if isempty(dicron),return;end897sigLen=length(dicron);898for diffIndex=1:sigLen899 if dicron(diffIndex)~=0900 tempp(diffIndex)=dicron(diffIndex)-od;901 else902 tempp(diffIndex)=0;903 end904end905clear dicron;906dicron=tempp;907clear tempp;908end909 910function [mini,minip,maxi,maxip]=seeklocales(tempsig,tempbegin,tempend)911tempMin=tempsig(tempbegin);912tempMax=tempsig(tempbegin);913minip=tempbegin;914maxip=tempbegin;915for j=tempbegin:tempend916 if tempsig(j)>tempMax917 tempMax=tempsig(j);918 maxip=j;919 elseif tempsig(j)<tempMin920 tempMin=tempsig(j);921 minip=j;922 end923end924 925mini=tempMin;926maxi=tempMax;927end928 929function [dicron]=seekdicrotic(tempdiff)930izcMin=0;931izcMax=0;932itemp=3;933tempLen=length(tempdiff)-3;934 935dicron=0;936 937tempdiff=smooth(tempdiff);938 939while itemp<=tempLen940 if (tempdiff(itemp)*tempdiff(itemp+1))<=0941 if tempdiff(itemp-2)<0942 if tempdiff(itemp+2)>=0943 izcMin=izcMin+1;944 tzcMin(izcMin)=itemp;945 end946 end947 948% if tempdiff(itemp-2)>0949% if tempdiff(itemp+2)<=0950% izcMax=izcMax+1;951% tzcMax(izcMax)=itemp;952% end953% end954 end955 956 itemp=itemp+1;957end958 959if izcMin==0 %big inflection960 itemp=3;961 tempMin=tempdiff(itemp);962 itempMin=itemp;963 964 while itemp<tempLen965 if tempdiff(itemp)<tempMin966 tempMin=tempdiff(itemp);967 itempMin=itemp;968 end969 itemp=itemp+1;970 end971 972 itemp=itempMin+1;973 while itemp<tempLen974 if tempdiff(itemp+1)<=tempdiff(itemp-1)975 dicron=itemp;976 return;977 end978 itemp=itemp+1;979 end980elseif izcMin==1981 dicron=tzcMin(izcMin);982 return;983else984 itemp=tzcMin(1);985 tempMax=tempdiff(itemp);986 itempMax=itemp;987 988 while itemp<tempLen989 if tempdiff(itemp)>tempMax990 tempMax=tempdiff(itemp);991 itempMax=itemp;992 end993 itemp=itemp+1;994 end995 996 for itemp=izcMin:-1:1997 if tzcMin(itemp)<itempMax998 dicron=tzcMin(itemp);999 return;1000 end1001 end1002end1003end1004 1005function [diap]=seekdiap(tempabp)1006diap=0;1007 1008[tt,ti]=max(tempabp);1009if ti==01010 diap=floor(length(tempabp)./2);1011else1012 diap=ti;1013end1014end1015","MATLAB"
1016"Biosensors","Cassey2016/PPG_Peak_Detection","func/my_func_standardizing_PPG.m",".m","422","11","function PPG_buffer = my_func_standardizing_PPG(PPG_buffer)1017% Standardizing PPG into zero-mean and uni-variance.1018 var_sig_PPG = var(PPG_buffer);1019 if var_sig_PPG == 01020 univar_sig_PPG = PPG_buffer;1021 else1022 univar_sig_PPG = sqrt(1/var_sig_PPG) * PPG_buffer;1023 end1024 zeromean_sig_PPG = univar_sig_PPG - mean(univar_sig_PPG);1025 PPG_buffer = zeromean_sig_PPG; % univariance for PPG 30 sec segment1026end","MATLAB"
1027"Biosensors","Cassey2016/PPG_Peak_Detection","func/my_func_prep_PPG_buffer.m",".m","554","17","function [PPG_buffer,fs_PPG] = my_func_prep_PPG_buffer(PPG_raw_buffer,fs_PPG)1028 % Resample PPG to 50 Hz.1029 if fs_PPG ~= 50 % Hz1030 PPG_down = resample(PPG_raw_buffer,50,fs_PPG);1031 fs_PPG = 50;1032 else1033 PPG_down = PPG_raw_buffer;1034 end1035 1036 PPG_buffer = PPG_down(:); % Make sure PPG is column vector1037 % Standardizing PPG in sub-function.1038 PPG_buffer = my_func_standardizing_PPG(PPG_buffer);1039 1040 % Filter signal.1041 [b, a] = butter(6,[0.5 20]/(fs_PPG/2)); % Bandpass filter.1042 PPG_buffer = filtfilt(b, a, PPG_buffer);1043end","MATLAB"
1044"Biosensors","Cassey2016/PPG_Peak_Detection","method_06/my_revise_run_wabp.m",".m","4330","109","function [r,ssf,my_avg0,A] = my_revise_run_wabp(abp,fs_abp)1045% Below was copied from Erick Andres Perez Alday's Github repository 1046% "" physionetchallenges / matlab-classifier-2020 "": 1047% https://github.com/physionetchallenges/matlab-classifier-2020/blob/master/Tools/PhysioNet-Cardiovascular-Signal-Toolbox-master/Tools/BP_Tools/run_wabp.m1048% WABP ABP waveform onset detector.1049% r = run_wabp(abp) obtains the onset time (in samples) 1050% of each beat in the ABP waveform.1051%1052% In: ABP (125Hz sampled)1053% Out: Onset sample time1054% 1055% Usage:1056% - ABP waveform must have units of mmHg1057%1058% Written by James Sun (xinsun@mit.edu) on Nov 19, 2005. This ABP onset1059% detector is adapted from Dr. Wei Zong's wabp.c.1060%1061% LICENSE: 1062% This software is offered freely and without warranty under 1063% the GNU (v3 or later) public license. See license file for1064% more information1065 1066% Dong changed: input should be 250 Hz for filtering.1067%% Input checks1068 % if nargin ~=11069 % error('exactly 1 argment needed');1070 % end1071 1072 if size(abp,2)~=11073 error('Input must be a <nx1> vector');1074 end1075 1076%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1077 1078 % scale physiologic ABP1079 offset = 1600;1080 scale = 20;1081 Araw = abp*scale-offset;1082 1083 % LPF1084 A = filter([1 0 0 0 0 -2 0 0 0 0 1],[1 -2 1],Araw)/24+30;1085 A = (A(4:end)+offset)/scale; % Takes care of 4 sample group delay1086 1087 % ------- Dong changed this: -------1088 A = A ./ std(A); % normalizing data is very important for my peak detection.1089 A = A - mean(A);1090 1091 % Slope-sum function1092 dypos = diff(A);1093 dypos(dypos<0) = 0;1094 % ssf = [0; 0; conv(ones(16,1),dypos)];1095 w = 16/125*fs_abp; % 125 Hz to 250 Hz.1096 ssf = [0; 0; conv(ones(w,1),dypos)];1097 1098 % Decision rule1099 first_8sec = 8*fs_abp;1100 % avg0 = sum(ssf(1:1000))/1000; % average of 1st 8 seconds (1000 samples) of SSF1101 avg0 = sum(ssf(1:first_8sec))/first_8sec;1102 Threshold0 = 3*avg0; % initial decision threshold1103 1104 % ignoring ""learning period"" for now1105 lockout = 0; % lockout >0 means we are in refractory1106 timer = 0;1107 % z = zeros(100000,1);1108 z = zeros(fs_abp*800,1);1109 counter = 0;1110 1111 % Dong: copied from wabp.c, 02/27/2020. % Dong change here. 02/27/2020.1112 TmDEF = 0.25; %5;% Dong change here. 02/27/2020.1113 max_min_thres = 0.1; %10;% Dong change here. 02/27/2020.1114 my_avg0 = zeros(size(abp));% Dong change here. 02/27/2020.1115 step_adjust_thres = 0.025; % it was 0.1 % Dong change here. 02/27/2020.1116 % for t = 50:length(ssf)-171117 for t = round(0.4*fs_abp):length(ssf)-w-11118 lockout = lockout - 1;1119 timer = timer + 1; % Timer used for counting time after previous ABP pulse1120 1121 if (lockout<1) & (ssf(t)>avg0+TmDEF) %(ssf(t)>avg0+5) % Not in refractory and SSF has exceeded threshold here % Dong change here. 02/27/2020.1122 timer = 0;1123 maxSSF = max(ssf(t:t+w)); % Find local max of SSF1124 minSSF = min(ssf(t-w:t)); % Find local min of SSF1125 if maxSSF > (minSSF + max_min_thres) %(minSSF + 10)% Dong change here. 02/27/2020.1126 onset = 0.01*maxSSF ; % Onset is at the time in which local SSF just exceeds 0.01*maxSSF1127 1128 tt = t-w:t;1129 dssf = ssf(tt) - ssf(tt-1);1130 BeatTime = find(dssf<onset,1,'last')+t-w-1;1131 counter = counter+1;1132 1133 if isempty(BeatTime)1134 counter = counter-1;1135 else1136 z(counter) = BeatTime;1137 end1138 Threshold0 = Threshold0 + step_adjust_thres*(maxSSF - Threshold0); % adjust threshold1139 avg0 = Threshold0 / 3; % adjust avg1140 1141 lockout = round(32/125*fs_abp); % lock so prevent sensing right after detection (refractory period)1142 end1143 end1144 1145 if timer > round(312/125*fs_abp) % Lower threshold if no pulse detection for a while1146 Threshold0 = Threshold0 - 0.1; %Threshold0 - 1; % Dong change here. 02/27/2020.1147 avg0 = Threshold0/3;1148 end1149 my_avg0(t,1) = avg0+TmDEF; % % Dong change here. 02/27/2020.1150 end1151 r = z(find(z))-2;1152end","MATLAB"
1153"Biosensors","Cassey2016/PPG_Peak_Detection","method_06/my_Elgendi_2013_method_III_peakdet.m",".m","1028","22","function [output_Elgendi_3_2013] = my_Elgendi_2013_method_III_peakdet(raw_PPG,fs_PPG)1154% -------------------------------------------------------------------------1155% This peak detection function was mentioned in this paper:1156% Elgendi, Mohamed, et al. 1157% ""Systolic peak detection in acceleration photoplethysmograms measured from 1158% emergency responders in tropical conditions."" PLoS One 8.10 (2013): e76585.1159% 1160 [r,ssf,my_avg0,A] = my_revise_run_wabp(raw_PPG,fs_PPG);1161% -------------------------------------------------------------------------1162 if isempty(r)1163 HR_Elgendi_3_2013 = 0; % there is no peak location.1164 r = 1;1165 else1166 HR_Elgendi_3_2013 = 60 * fs_PPG ./ diff(r); % calculate the HR.1167 end1168 A = [A;0;0;0;]; % add zero1169 A(1:6) = A(7); % first six plots are all high amplitude.1170 output_Elgendi_3_2013 = struct('PPG_peak_loc_Elgendi_3_2013',r,...1171 'HR_Elgendi_3_2013',HR_Elgendi_3_2013,...1172 'filtered_PPG_Elgendi_3_2013',A,...1173 'thres_Elgendi_3_2013',my_avg0);1174end","MATLAB"
1175"Biosensors","Cassey2016/PPG_Peak_Detection","method_03_and_04/my_Elgendi_2013_method_I_peakdet.m",".m","3762","111","function [output_Elgendi_1_2013] = my_Elgendi_2013_method_I_peakdet(raw_PPG, delta, fs_PPG)1176% -------------------------------------------------------------------------1177% Dong add this on 02/25/2020, based on this paper:1178% Elgendi, Mohamed, et al. 1179% ""Systolic peak detection in acceleration photoplethysmograms measured from 1180% emergency responders in tropical conditions."" PLoS One 8.10 (2013): e76585.1181% 1182% (1): bandpass filter (0.5-8Hz)1183 [b, a] = butter(6,[0.5 8]/(fs_PPG/2)); % bandpass filter 0.5-10Hz, changed from 0.5-20 to 0.5-9 Hz at 11/21/20181184 raw_PPG = filtfilt(b, a, raw_PPG); % -> AC component1185 raw_PPG = raw_PPG ./ std(raw_PPG); % normalizing data is very important for my peak detection.1186 raw_PPG = raw_PPG - mean(raw_PPG);1187 1188 debugging_plot_flag = false; % only for plotting debugging figures. 1189% -------------------------------------------------------------------------1190% Below code is copied from: http://billauer.co.il/peakdet.html1191% PEAKDET Detect peaks in a vector1192% [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local1193% maxima and minima (""peaks"") in the vector V.1194% MAXTAB and MINTAB consists of two columns. Column 11195% contains indices in V, and column 2 the found values.1196% 1197% With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices1198% in MAXTAB and MINTAB are replaced with the corresponding1199% X-values.1200%