Skip to content

Commit

Permalink
new large spike compensation
Browse files Browse the repository at this point in the history 
git-svn-id: https://svn.code.sf.net/p/pspm/svn/trunk@527 23ec9670-3a1e-4a12-9187-a0ba8bce340d
  • Loading branch information
Tobias Moser committed Feb 16, 2018
1 parent 65d90ee commit 25caba1
Showing 1 changed file with 84 additions and 51 deletions.
135 changes: 84 additions & 51 deletions pspm_ppu2hb.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -25,91 +29,145 @@
sts = -1;
outinfo = struct();
global settings;
if isempty(settings), pspm_init; end;
if isempty(settings), pspm_init; end

%% check input
% -------------------------------------------------------------------------
if nargin < 1
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
% -------------------------------------------------------------------------
fprintf('Heartbeat detection for %s ... \n', fn);

% 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
if length(t_ppu) ~= length(ppu_corr)
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(:);
Expand All @@ -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
Expand All @@ -139,29 +198,3 @@
return;

end


























0 comments on commit 25caba1

Please sign in to comment.