function [EEG, PERCBEATSOK, SDIBI, MEDIANIBI] = mig_ContECGToIBI(EEG, ECGchan, lowerbeatthresh, upperbeatthresh, useIBIptsOrBPM, checkitout, createnewevents) %mig_ContECGtoIBI: converts raw ECG channel into continuous IBI data % note: this bit of code requires the excellent "peakfinder" algorithm (on % matlab central) % added: initial peakfinder now uses an area of 50 ms to define what counts % as a peak (mainly so that sudden problems with noisy data (falling % electrodes..) will not dominate the peakfinder. % % IN: EEG, ECGchan, lowerbeatthresh, upperbeatthresh, useIBIptsOrBPM, checkitout, createnewevents % OUT: EEG, PERCBEATSOK(percentage ok IBI), SDIBI (std of IBI), MEDIANIBI % (median of IBI) % % inputs: % EEG: EEG struct in % ECGchan: which channel concerns ECG? % outputs: % lower beat threshold. Threshold for any type of beat. If function starts % to go wrong, use a higher lower. Otherwise enter 0! Normally 50% of avg % peak. % upper beat threshold. Threshold that should sit somewhere between the % "up" and "down" beat of the ECG. Can be tweaked for people who have % tricky heartrates (e.g. with big "down" beat) % useIBIptsOrBPM: set to 1 if output needs to be beats per minute as % opposed to (default) IBI (ms) % checkitout: set to 1 if a quick check is desired. Shows 1) after filtering, % 2) initial peak detection, 3) after continuous conversion. % createnewevents: set to 1 if markers for beats are needed (takes a % while!) % outputs % EEG struct. % % example: EEG = mig_ContECGToIBI(pop_loadset('mydataset.set'), 4, 0, 0, 1, 1, 0) % % Copyright M. Spapé, 2015. % 7/4/2015: Fix for incorrect treatment of flipped data. % 7/4/2015: Fixed IBI / BPM conversion messing up EEG struct. % 7/4/2015: Changed default criterion for peak threshold. Now at 0.75 * the 75th percentile of peaks. Less sensitive to extremely bad data. % 20/8/2015: Added initial check for really horrible data. Algorithm goes % through every 5 s of data and entirely removes the data if the SD of data % is larger than median+20*inter quartile range. threshfortinyIBI = 30 / (1000 / EEG.srate) ; %we imagine the same beat is detected twice if the IBI is 30 reallyhighIBI = 1500 / (1000 / EEG.srate); reallyhighSDIBI = 150 / (1000 / EEG.srate); reallylowIBI = 400 / (1000 / EEG.srate); suddenbigjumpinIBI = 300 / (1000 / EEG.srate); initpeakfinderarea = 100 / (1000 / EEG.srate); DATAWAS = EEG; EEG = pop_select(EEG, 'channel', ECGchan); if std(EEG.data) == 0, PERCBEATSOK = 0; SDIBI = 10000; MEDIANIBI = 0; else EEG = pop_eegfiltnew(EEG, 2, 0); %takes about 4 minutes :( EEG = pop_eegfiltnew(EEG, 0, 100); %filter everything out above 40 hz. EEG = pop_eegfiltnew(EEG, 46, 54, [], 1, [], 0);%filter everything out between 46 and 54; tmp_allsds = []; for ctime = 2501:5000:size(EEG.data,2)-2500, tmp_allsds = [tmp_allsds std(EEG.data(1,ctime-2500:ctime+2499))]; end horrordatathreshold = median(tmp_allsds)+20*(prctile(tmp_allsds,75)-prctile(tmp_allsds,25)); tmp_removedhowmuch = 0; for ctime = 2501:5000:size(EEG.data,2)-2500, if std(EEG.data(1,ctime-2500:ctime+2499))>horrordatathreshold, EEG.data(1,ctime-2500:ctime+2499) = 0; tmp_removedhowmuch = tmp_removedhowmuch + 5; end end disp (['removed ' num2str(tmp_removedhowmuch) 's of data from view.']); if checkitout, pop_eegplot(EEG, 1, 0,0); end; disp ('Finding positive peaks.'); [~, peakmagsplus] = peakfinder(EEG.data, initpeakfinderarea); peakmagsplus(peakmagsplus<0) = []; disp ('Finding negative peaks.'); [~, peakmagsmin] = peakfinder(-EEG.data, initpeakfinderarea); peakmagsmin(peakmagsmin<0) = []; if prctile(peakmagsmin,75)>prctile(peakmagsplus,75), disp ('Higher negative than positive peak magnitudes. Flipping data!'); EEG.data = -EEG.data; peakmags = peakmagsmin; else disp ('Higher postive than negative peak magnitudes. Keeping data!'); peakmags = peakmagsplus; end if upperbeatthresh == 0, upperbeatthresh = 0.75*prctile(peakmags,75); disp (['Setting upper beat threshold to 75% of 75 percentile of peak mag: ' num2str(upperbeatthresh) '.']); end allpeaks = zeros(1,3); for ctime = 1:512:(size(EEG.data,2)-3999) [peaklocs, peakmags] = peakfinder(EEG.data(ctime:ctime+3999), [], upperbeatthresh); allpeaks = [allpeaks; (peaklocs + (ctime-1))' peakmags' (zeros(size(peaklocs)) + std(EEG.data(ctime:ctime+3999)))']; end allpeaks(1,:) = []; [~, keepthese] = unique(allpeaks(:,1)); allpeaks = allpeaks(keepthese,:); clear keepthese; allpeaks = [allpeaks [allpeaks(1,1); diff(allpeaks(:,1))]]; allpeaks = [allpeaks ones(size(allpeaks,1),1)]; %get rid of really noisy data! allpeaks(1:length(allpeaks),[5 6])=1; allpeaks(unique([find(allpeaks(:,3) > 2* median(allpeaks(:,3))); find(allpeaks(:,3) > 2* median(allpeaks(:,3)))-1;find(allpeaks(:,3) > 2* median(allpeaks(:,3)))+1]),6) = 0; if lowerbeatthresh == 0, threshforlowerbeat = median(allpeaks(find(allpeaks(allpeaks(:,6)==1)),2))/2; %define weird peak magnitude as those with > .5 * median else threshforlowerbeat = lowerbeatthresh; end %we kick out everything with low magnitude and extrapolate IBI to %high magnitudes for dotimes = 1:3 for cpeak = 2:size(allpeaks,1) if allpeaks(cpeak,2) < threshforlowerbeat, if allpeaks(cpeak,2) > allpeaks(cpeak-1,2), allpeaks(cpeak-1,5) = 0; else allpeaks(cpeak,5) = 0; end end end allpeaks = allpeaks(allpeaks(:,5) == 1,:); allpeaks(:,4) = [allpeaks(1,1); diff(allpeaks(:,1))]; threshforlowerbeat = median(allpeaks(find(allpeaks(allpeaks(:,6)==1)),2))/2; end ; for cpeak = 2:size(allpeaks,1) if allpeaks(cpeak,4) < threshfortinyIBI, if allpeaks(cpeak,2) > allpeaks(cpeak-1,2), allpeaks(cpeak-1,5) = 0; else allpeaks(cpeak,5) = 0; end end end allpeaks = allpeaks(allpeaks(:,5) == 1,:); allpeaks(:,4) = [allpeaks(1,1); diff(allpeaks(:,1))]; SDIBI = std(diff(allpeaks(allpeaks(:,4) reallyhighSDIBI, for dotimes = 1:3 disp(['WARNING! SD delta(IBI) found rather high (' num2str(SDIBI) '). Countering with detection based on IBI']); IBITHRESH = median(allpeaks(find(allpeaks(:,4) allpeaks(cpeak-1,2), allpeaks(cpeak-1,5) = 0; else allpeaks(cpeak,5) = 0; end end end allpeaks = allpeaks(allpeaks(:,5) == 1,:); allpeaks(:,4) = [allpeaks(1,1); diff(allpeaks(:,1))]; IBITHRESH = median(allpeaks(find(allpeaks(:,4)reallyhighIBI),5) = 0; allpeaks(find(allpeaks(:,4)suddenbigjumpinIBI)+1,7) = 0; allpeaks(:,8) = allpeaks(:,5) .* allpeaks(:,6) .* allpeaks(:,7); TOOHIGH = prctile(allpeaks(find(allpeaks(:,8)==1),4),50) + 2* (prctile(allpeaks(find(allpeaks(:,8)==1),4),75) - prctile(allpeaks(find(allpeaks(:,8)==1),4),25)); TOOLOW = prctile(allpeaks(find(allpeaks(:,8)==1),4),50) - 2* (prctile(allpeaks(find(allpeaks(:,8)==1),4),75) - prctile(allpeaks(find(allpeaks(:,8)==1),4),25)); allpeaks(find(allpeaks(:,4)>TOOHIGH),5) = 0; allpeaks(find(allpeaks(:,4)