From 25caba15f86d888ece9b5a6ab87075c45927f4f9 Mon Sep 17 00:00:00 2001 From: Tobias Moser Date: Fri, 16 Feb 2018 12:49:45 +0000 Subject: [PATCH] new large spike compensation git-svn-id: https://svn.code.sf.net/p/pspm/svn/trunk@527 23ec9670-3a1e-4a12-9187-a0ba8bce340d --- pspm_ppu2hb.m | 135 +++++++++++++++++++++++++++++++------------------- 1 file changed, 84 insertions(+), 51 deletions(-) diff --git a/pspm_ppu2hb.m b/pspm_ppu2hb.m index 0497139ca..5ec31ca92 100644 --- a/pspm_ppu2hb.m +++ b/pspm_ppu2hb.m @@ -12,6 +12,10 @@ % diagnostics : [true/FALSE] displays some debugging information % replace : [true/FALSE] replace existing heartbeat channel. % If multiple channels are present, replaces last. +% lsm : [integer] large spikes mode compensates for +% large spikes while generating template by +% removing the [integer] largest percentile of +% spikes from consideration %__________________________________________________________________________ % PsPM 3.1 % (C) 2016 Samuel Gerster (University of Zurich), Tobias Moser (University of Zurich) @@ -25,7 +29,7 @@ sts = -1; outinfo = struct(); global settings; -if isempty(settings), pspm_init; end; +if isempty(settings), pspm_init; end %% check input % ------------------------------------------------------------------------- @@ -33,19 +37,21 @@ warning('ID:invalid_input', 'No input. Don''t know what to do.'); return; elseif ~ischar(fn) warning('ID:invalid_input', 'Need file name string as first input.'); return; -elseif nargin < 2 +elseif nargin < 2 || isempty(chan) chan = 'ppu'; elseif ~isnumeric(chan) && ~strcmp(chan,'ppu') warning('ID:invalid_input', 'Channel number must be numeric'); return; -end; +end %%% Process options % Display diagnostic plots? default is "false" -try if ~islogical(options.diagnostics),options.diagnostics = false;end; - catch, options.diagnostics = false; end; +try if ~islogical(options.diagnostics),options.diagnostics = false;end + catch, options.diagnostics = false; end % Replace existing heartbeat channel? default is "false" -try if ~islogical(options.replace),options.replace = false;end; - catch, options.replace = false; end; +try if ~islogical(options.replace),options.replace = false;end + catch, options.replace = false; end +try if ~isnumeric(options.lsm),options.lsm = 0;end + catch, options.lsm = 0; end %% user output % ------------------------------------------------------------------------- @@ -53,48 +59,92 @@ % get data % ------------------------------------------------------------------------- -[nsts, infos, data] = pspm_load_data(fn, chan); -if nsts == -1, return; end; +[nsts, ~, data] = pspm_load_data(fn, chan); +if nsts == -1, return; end if numel(data) > 1 fprintf('There is more than one PPU channel in the data file. Only the first of these will be analysed.'); data = data(1); -end; - -if ~any(strcmp(data{1,1}.header.chantype,{'ppu'})) +end +% Check that channel is ppu +if ~strcmp(data{1,1}.header.chantype,'ppu') warning('ID:not_allowed_channeltype', 'Specified channel is not a PPU channel. Don''t know what to do!') return; end +%% Large spikes mode +%-------------------------------------------------------------------------- +ppu = data{1}.data; +% large spike mode +if options.lsm + fprintf('Entering large spikes mode. This might take some time.'); + % Look for all peaks lower than 200 bpm (multiple of two in heart rate + % to compensate for absolute value and therefore twice as mani maxima) + [pks,pis] = findpeaks(abs(ppu),... + 'MinPeakDistance',30/200*data{1}.header.sr); + % Ensure at least one spike is removed by adapting quantil to realistic + % values, given number of detected spikes + q = floor(length(pks)*(1-options.lsm/100))/length(pks); + % define large spikes index as last lsm percentile (or as adapted above) + lsi = pks>quantile(pks,q); + %define a minimum peak prominence 2/3 of non large spikes range (more + %or less) + minProm = max(pks(~lsi))*2/3; + % save indexes of large spikes for later removal while generating + % template + lsi = pis(lsi); + fprintf(' done.\n'); +else + minProm = range(ppu)/3; +end + %% Create template %-------------------------------------------------------------------------- +fprintf('Creating template. This might take some time.'); % Find prominent peaks for a max heart rate of 200 bpm -[~,pi] = findpeaks(data{1}.data,... +[~,pis] = findpeaks(data{1}.data,... 'MinPeakDistance',60/200*data{1}.header.sr,... - 'MinPeakProminence',range(data{1}.data/3)); -if isempty(pi) - warning('ID:NoPulse', 'No pulse found, nothing done.'); return; + 'MinPeakProminence',minProm); + +if options.lsm + % Remove large spikes from + [~,lsi_in_pis,~] = intersect(pis,lsi); + pis(lsi_in_pis) = []; end -% get pulse epochs lower limit -d = min(diff(pi)); -epochs_l = floor(pi(2:end-1)-.3*d); +% handle possible errors +if isempty(pis),warning('ID:NoPulse', 'No pulse found, nothing done.');return;end + +% get pulse period lower limit (assumed onset) as 30% of smalest period +% before detected peaks +min_pulse_period = min(diff(pis)); +period_index_lower_bound = floor(pis(2:end-1)-.3*min_pulse_period); +fprintf('...'); -%Create template from mean of peak time-locked ppu pulse epochs -pulses = cell2mat(arrayfun(@(x) data{1}.data(x:x+d),epochs_l','un',0)); +% Create template from mean of peak time-locked ppu pulse periods +pulses = cell2mat(arrayfun(@(x) data{1}.data(x:x+min_pulse_period),period_index_lower_bound','un',0)); template = mean(pulses,2); +fprintf('done.\n'); + +% handle diagnostic plots relevant to template building if options.diagnostics t_template = (0:length(template)-1)'/data{1}.header.sr; - t_pulses = repmat(t_template,1,length(pi)-2); + t_pulses = repmat(t_template,1,length(pis)-2); figure - plot(t_pulses,pulses,'--',t_template,template,'k','lineWidth',2) + plot(t_pulses,pulses,'--') + set(gca,'NextPlot','add') + plot(t_template,template,'k','lineWidth',3) + xlabel('time [s]') + ylabel('Amplitude') + title('Generated ppu template (bk) and pulses used (colored)') end %% Cross correlate the signal with the template and find peaks %-------------------------------------------------------------------------- +fprintf('Applying template.'); ppu_corr = xcorr(data{1}.data,template)/sum(template); % Truncate ppu_xcorr and realigne it so the max correlation corresponds to % templates peak and not beginning of template. -ppu_corr = ppu_corr(length(data{1}.data)-floor(.3*d):end-floor(.3*d)); +ppu_corr = ppu_corr(length(data{1}.data)-floor(.3*min_pulse_period):end-floor(.3*min_pulse_period)); if options.diagnostics t_ppu = (0:length(data{1}.data)-1)'/data{1}.header.sr; figure @@ -102,14 +152,22 @@ length(t_ppu) end plot(t_ppu,ppu_corr,t_ppu,data{1}.data) + xlabel('time [s]') + ylabel('Amplitude') + title('ppu cross-corelated with template and ppu') + legend('ppu (X) template','ppu') end % Get peaks that are at least one template width appart. These are the best % correlation points. -[~,hb] = findpeaks(ppu_corr/max(ppu_corr),data{1}.header.sr,'MinPeakdistance',d/data{1}.header.sr); +[~,hb] = findpeaks(ppu_corr/max(ppu_corr),... + data{1}.header.sr,... + 'MinPeakdistance',min_pulse_period/data{1}.header.sr); +fprintf(' done.\n'); %% Prepare output and save %-------------------------------------------------------------------------- % save data +fprintf('Saving data.'); msg = sprintf('Heart beat detection from ppu with cross correlation HB-timeseries added to data on %s', date); newdata.data = hb(:); @@ -120,6 +178,7 @@ write_options = struct(); write_options.msg = msg; +% Replace last existing channel or save as new channel if options.replace write_action = 'replace'; else @@ -139,29 +198,3 @@ return; end - - - - - - - - - - - - - - - - - - - - - - - - - -