function PulseComplete(name,type,TRi,HRrate,varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Pulse Covariate generator % % Inputs: % type - Either 'single', 'multi', or 'both' % This term specifies whether a single concatenated output or % multiple outputs are desired % TR - the TR interval used in for the scan % HRrate - 0 for high average heartrate (>100), 1 for low average HR rate (<100) % varargin - series of file names of the format 'xxxx.ref' % Note: if a file is missing, a place holder of ' ' MUST BE USED % % Derived from heavily modified PhLEM toolbox code % Written by Omar H Butt Aug 2010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if HRrate == 0 HRrate = 'high'; elseif HRrate == 1 HRrate = 'low'; end HRrate = lower(HRrate); % if mod(n,20)~=0; % disp ('ERROR: Selected shift must be divisible by 20ms'); % else n = 0; % Hard coded shift desired. Note, must be in multiples of 20ms. spacer = n / 20; CO = []; gotone=0; namer = 1; totallength=0; for i = 2:2:length(varargin) totallength = totallength + varargin{i}; end HRs =[]; LRs =[]; COS =[]; baseset =[]; numb = 1; switch lower(type) case 'single' for i = 1:2:(length(varargin)-1) isempty = strcmp(varargin{i},' '); if isempty == 1 currentoutput = zeros(8,varargin{i+1}); else if gotone == 0 namer = i; gotone = 1; end currentfile = varargin{i}; TR = varargin{i+1}; [ currentoutput HR LR ] = PulseSplitCore(currentfile,TR,spacer,TRi,HRrate); HRs = [HRs HR]; LRs = [LRs LR]; end COS = [COS currentoutput]; end COS = COS'; fname = name; HRF = mean(HRs); LRF = mean(LRs); writeme(COS,n,fname,HRF,LRF); disp('Complete!') case 'multi' HRx = 0; LRx = 0; for i = 1:2:(length(varargin)-1) isempty = strcmp(varargin{i},' '); if isempty == 1 currentoutput = zeros(8,varargin{i+1}); else if gotone == 0 namer = i; gotone = 1; end currentfile = varargin{i}; TR = varargin{i+1}; [ currentoutput HRx LRx ] = PulseSplitCore(currentfile,TR,spacer,TRi,HRrate); end if i == 1 desiredout = currentoutput; else desiredout = [baseset currentoutput]; end desiredout(1,end+1:totallength)=0; fname = [sprintf('Scan%u-',numb) name]; writeme(desiredout',n,fname,HRx,LRx); baseset = [baseset zeros(8,varargin{i+1})]; numb = numb+1; end disp('Complete!') case 'both' baseset = []; HRx = 0; LRx = 0; for i = 1:2:(length(varargin)-1) isempty = strcmp(varargin{i},' '); if isempty == 1 currentoutput = zeros(8,varargin{i+1}); else if gotone == 0 namer = i; gotone = 1; end currentfile = varargin{i}; TR = varargin{i+1}; [ currentoutput HRx LRx ] = PulseSplitCore(currentfile,TR,spacer,TRi,HRrate); end if i == 1 desiredout = currentoutput; else desiredout = [baseset currentoutput]; end desiredout(1,end+1:totallength)=0; fname = [sprintf('Scan%u-',numb) name]; writeme(desiredout',n,fname,HRx,LRx); baseset = [baseset zeros(8,varargin{i+1})]; COS = [COS currentoutput]; numb = numb+1; end COS = COS'; fname = name; HRF = mean(HRs); LRF = mean(LRs); writeme(COS,n,fname,HRF,LRF); disp('Complete!') otherwise disp('Unknown desired output parameters.') end % end function writeme (CO,n,fname,HRF,LRF) Hstring = [ ';VB98'; ';REF1'; '; '; ]; ORs = sprintf('; Cardiac Rate: %f',HRF); Ms = sprintf('; File length: %u',length(CO)); writetotext(CO(:,1),[sprintf('H-Sin1-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,2),[sprintf('H-Cos1-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,3),[sprintf('H-Sin2-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,4),[sprintf('H-Cos2-(%u)-',n) fname],Hstring,ORs,Ms) ORs = sprintf('; Resp Rate: %f',LRF); writetotext(CO(:,5),[sprintf('L-Sin1-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,6),[sprintf('L-Cos1-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,7),[sprintf('L-Sin2-(%u)-',n) fname],Hstring,ORs,Ms) writetotext(CO(:,8),[sprintf('L-Cos2-(%u)-',n) fname],Hstring,ORs,Ms) function [ genoutput Highrate Lowrate ] = PulseSplitCore(fname,TRtotal,spacer,TR,HRrate) Hz = 50; TR = TR*1000; %%%Text Analysis Portion below. Adapted for .ref file manipulation fclose('all');fid=fopen(fname); linestoignore = 9; ignore=textscan(fid,'%s',linestoignore,'delimiter','\n'); %Ignore first 2 rows. [ data ] = textscan(fid,'%n'); %Read data until end of data. signal = data{1}; fclose(fid); signal=signal'; % Find Proper length of total pulse array, and interpolate it down to the % new shorter size to accomodate for the slightly longer TR interval sampR=1000/Hz; arraylength=TRtotal*TR/sampR; PSig = sizecorrection(signal,arraylength); signal = PSig; Time = (1:length(signal))./Hz; %%%Apply smoothing/filtering to the signal [HIGHsignal Hevents] = PhLEM_H(signal,Hz,HRrate); [LOWsignal Levents] = PhLEM_L(signal,Hz,HRrate); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Highes = PulseRegress(Hevents,Time,TR,spacer); Lowes = PulseRegress(Levents,Time,TR,spacer); Highrate = length(find(Hevents)) / Time(end) * 60; disp([fname sprintf(' Cardiac Rate: %f',Highrate)]) Lowrate = length(find(Levents)) / Time(end) * 60; disp([fname sprintf(' Resp Rate: %f',Lowrate)]) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Output genoutput = [Highes Lowes]; genoutput=genoutput'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Resize function to make correct length based on TR multiple %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [outputsignal] = sizecorrection(signal,newlength) x = 1:length(signal); shifter = length(signal)/newlength; outputsignal = zeros (1,newlength); currentspot = 1; for cc=1:newlength if currentspot > length(signal) currentspot = length(signal); else outputsignal(1,cc) = interp1(x,signal,currentspot,'linear'); currentspot = currentspot + shifter; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % delta = search rate for looking for signal peaks % peak_rise = amplitude threshold for looking for peaks % filter = finter or smoothing type ('butter' or 'gaussian') % Wn = filter/smoothing kernel. If bandpass filter, % normalized frequencies. If gaussian smoothing, % then kernel fwhm. % Hz = Frequency of the sampled data signal %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Filter and smoothing function for HIGH frequencies %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x events] = PhLEM_H(TS,Hz, HRrate) if strcmp(HRrate,'low') delta = Hz/3.5; Wn = [1 10]/Hz; else delta = Hz/6; % Default from Hz/6; Wn = [1 14]/Hz; % Defauly [ 1 12]; a range from 9 to 12 works for second number end peak_rise = 0.1; %Transform via ABS x = abs(TS); x = x - mean(x); % Use butter filter if ~isempty(Wn) [b a] = butter(1,Wn); xbutter = filter(b,a,x); end; % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % $$$ % -- Schlerf (Jul 24, 2008) % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Take a first pass at peak detection: mxtb = peakdet(x,1e-14); % set the threshold based on the 20th highest rather than the highest: sorted_peaks = sort(mxtb(:,2)); if length(sorted_peaks) > 20 peak_resp=sorted_peaks(end-20); else if length(sorted_peaks) > 15 peak_resp=sorted_peaks(end-15); else if length(sorted_peaks) > 10 peak_resp=sorted_peaks(end-10); else if length(sorted_peaks) > 5 peak_resp=sorted_peaks(end-5); else peak_resp = sorted_peaks(1); end end end end % And now do a second pass, more robustly filtered, to get the actual peaks: mxtb = peakdet(x,peak_rise*peak_resp); events = zeros(size(TS)); peaks = mxtb(:,1); dpeaks = diff(peaks); kppeaks = find(dpeaks > delta); newpeaks = peaks([1 kppeaks'+1]); events(newpeaks) = 1; % DEBUG CATCH: if length(newpeaks) < 5; warning('Program Crash: High Freq peaks hard to detect; try different HRrate field?') keyboard; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Filter and smoothing function for LOW frequencies %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x events] = PhLEM_L(TS,Hz,HRrate) if strcmp(HRrate,'low') delta = Hz*1.5; Wn = Hz*0.3; else delta = Hz*1.5; %Default is Hz*1.5 Wn = Hz*0.35; %Default is Hz*0.4 end peak_rise = 0.5; % Transform via ABS x = abs(TS); x = x - mean(x); % Use Gaussian Filter x = smooth_kernel(x(:),Wn); x = x'; % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % $$$ % -- Schlerf (Jul 24, 2008) % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Take a first pass at peak detection: [mxtb mntb] = peakdet(x,1e-14); % set the threshold based on the 20th highest rather than the highest: sorted_peaks = sort(mxtb(:,2)); if length(sorted_peaks) > 20 peak_resp=sorted_peaks(end-20); else if length(sorted_peaks) > 15 peak_resp=sorted_peaks(end-15); else if length(sorted_peaks) > 10 peak_resp=sorted_peaks(end-10); else if length(sorted_peaks) > 5 peak_resp=sorted_peaks(end-5); else peak_resp = sorted_peaks(1); end end end end % And now do a second pass, more robustly filtered, to get the actual peaks: [mxtb mntb] = peakdet(x,peak_rise*peak_resp); events = zeros(size(TS)); peaks = mxtb(:,1); dpeaks = diff(peaks); kppeaks = find(dpeaks > delta); newpeaks = peaks([1 kppeaks'+1]); events(newpeaks) = 1; % DEBUG CATCH: if length(newpeaks) < 5; delta = Hz*1.275; [mxtb mntb] = peakdet(x,1e-14); % set the threshold based on the 20th highest rather than the highest: sorted_peaks = sort(mxtb(:,2)); if length(sorted_peaks) > 20 peak_resp=sorted_peaks(end-20); else if length(sorted_peaks) > 15 peak_resp=sorted_peaks(end-15); else if length(sorted_peaks) > 10 peak_resp=sorted_peaks(end-10); else if length(sorted_peaks) > 5 peak_resp=sorted_peaks(end-5); else peak_resp = sorted_peaks(1); end end end end % And now do a second pass, more robustly filtered, to get the actual peaks: [mxtb mntb] = peakdet(x,peak_rise*peak_resp); events = zeros(size(TS)); peaks = mxtb(:,1); dpeaks = diff(peaks); kppeaks = find(dpeaks > delta); newpeaks = peaks([1 kppeaks'+1]); events(newpeaks) = 1; disp('Low Freq peaks hard to detect, switching delta') if length(newpeaks) < 5; warning('Program Crash: Low Freq peaks still hard to detect; Bad Data?') end % keyboard; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%To run a gaussian smoothing kernel over the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y=smooth_kernel(y,sigma) N=round(sigma*5)*2; x=[1:N]; v=normpdf(x,(N+1)/2,sigma); pb(1:length(v)-1)=y(1); pe(1:length(v)-1)=y(end); y=[pb';y;pe']; y=conv(y,v); cut=round(1.5*N); y=y(cut-1:end-cut+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Write to text file function %matx = matrix to output %N = number of sectors/rows %file = file name to output to %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function writetotext(M,file,Hstring,ORs,l) matx = (M); fid = fopen(file, 'wt'); % Open for writing for i=1:size(Hstring,1) fprintf(fid, '%s ', Hstring(i,:)); fprintf(fid, '\n'); end fprintf(fid, '%s ', ORs); fprintf(fid, '\n'); fprintf(fid, '%s ', l); fprintf(fid, '\n'); for n=1:5 fprintf(fid, ';'); fprintf(fid, '\n'); end for i=1:size(matx,1) fprintf(fid, '%d ', matx(i,:)); fprintf(fid, '\n'); end fclose(fid); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [maxtab, mintab]=peakdet(v, delta) %PEAKDET Detect peaks in a vector % [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local % maxima and minima ("peaks") in the vector V. % A point is considered a maximum peak if it has the maximal % value, and was preceded (to the left) by a value lower by % DELTA. MAXTAB and MINTAB consists of two columns. Column 1 % contains indices in V, and column 2 the found values. % Eli Billauer, 3.4.05 (Explicitly not copyrighted). % This function is released to the public domain; Any use is allowed. maxtab = []; mintab = []; v = v(:); % Just in case this wasn't a proper vector if (length(delta(:)))>1 error('Input argument DELTA must be a scalar'); end if delta <= 0 error('Input argument DELTA must be positive'); end mn = Inf; mx = -Inf; mnpos = NaN; mxpos = NaN; lookformax = 1; for i=1:length(v) this = v(i); if this > mx, mx = this; mxpos = i; end if this < mn, mn = this; mnpos = i; end if lookformax if this < mx-delta maxtab = [maxtab ; mxpos mx]; mn = this; mnpos = i; lookformax = 0; end else if this > mn+delta mintab = [mintab ; mnpos mn]; mx = this; mxpos = i; lookformax = 1; end end end function phase = Ph_expand(event_series,time) % function phase = PhLEM_phase_expand(event_series,time); % % Takes a series of events defined as 1's or 0's, and an index vector % for real time (in seconds), and returns the unwarped phase between % the events. Works by assuming that the distance between consecutive % events is 2pi. % % Inputs: % event_series = 1xN array of 1's (events) or 0's (non-events); % time = 1xN array of timestamps in real time (seconds); % % Written by T. Verstynen (November 2007) % % Liscensed under the GPL public license (version 3.0) % Time Stamps of Events events = find(event_series); % Phase Interpolation uwp_phase = [0:length(events)-1]; uwp_phase = 2*pi.*uwp_phase; phase = interp1(time(events),uwp_phase,time,'spline','extrap'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function [C] = PulseRegress(events,Time,TR,spacer) phase = [ 1 2 ]; for p = phase tmp_phase = Ph_expand(events,Time); eval(sprintf('dphase%d = %d*tmp_phase;',p,phase(p))); end; A = dphase1; B = dphase2; TR=TR/1000; % Setup TR markers bottom = mean(diff(Time)); ender= Time(end)-(TR/(2)); tr_timestamps = (0:TR:ender); tr_timestamps = tr_timestamps./bottom; %Analysis if (spacer ~= 0) && (spacer > 0) fillerA = A(1,end-spacer+1:end); fillerB = B(1,end-spacer+1:end); phasev1 = [ fillerA A(1,1:end-spacer) ]; phasev2 = [ fillerB B(1,1:end-spacer) ]; else if (spacer ~= 0) && (spacer < 0) fillerA = A(1,1:abs(spacer)); fillerB = B(1,1:abs(spacer)); phasev1 = [ A(1,1+abs(spacer):end) fillerA ]; phasev2 = [ B(1,1+abs(spacer):end) fillerB ]; else phasev1 = A; phasev2 = B; end end C = []; for exp = phase zz= eval(sprintf('phasev%d', phase(exp))); nphs = interp1(zz,tr_timestamps,'linear','extrap'); C = [C [sin(nphs)]']; C = [C [cos(nphs)]']; end