对脉搏波的波峰波谷以及重搏波拐点的检测

function [onsetp,peakp,dicron] = delineator(abpsig,abpfreq) % This program is intended to delineate the fiducial points of pulse waveforms % Inputs: % abpsig: input as original pulse wave signals; % abpfreq: input as the sampling frequency; % Outputs: % onsetp: output fiducial points as the beginning of each beat; % peakp: output fiducial points as systolic peaks; % dicron: output fiducial points as dicrotic notches; % Its delineation is based on the self-adaptation in pulse waveforms, but % not in the differentials. % Reference: % BN Li, MC Dong & MI Vai (2010) % On an automatic delineator for arterial blood pressure waveforms % Biomedical Signal Processing and Control 5(1) 76-81. % LI Bing Nan @ University of Macau, Feb 2007 % Revision 2.0.5, Apr 2009 %Initialization peakIndex=0; onsetIndex=0; dicroIndex=0; stepWin=2*abpfreq; closeWin=floor(0.1*abpfreq); %invalid for pulse beat > 200BPM sigLen=length(abpsig); peakp=[]; onsetp=[]; dicron=[]; %lowpass filter at first coh=25; %cutoff frequency is 25Hz coh=coh*2/abpfreq; od=3; %3rd order bessel filter [B,A]=besself(od,coh); abpsig=filter(B,A,abpsig); abpsig=10*abpsig; abpsig=smooth(abpsig); %Compute differentials ttp=diff(abpsig); diff1(2:sigLen)=ttp; diff1(1)=diff1(2); diff1=100*diff1; clear ttp; diff1=smooth(diff1); if sigLen>12*abpfreq tk=10; elseif sigLen>7*abpfreq tk=5; elseif sigLen>4*abpfreq tk=2; else tk=1; end %Seek average threshold in original signal if tk>1 %self-learning threshold with interval sampling tatom=floor(sigLen/(tk+2)); for ji=1:tk %search the slopes of abp waveforms sigIndex=ji*tatom; tempIndex=sigIndex+abpfreq; [tempMin,jk,tempMax,jl]=seeklocales(abpsig,sigIndex,tempIndex); tempTH(ji)=tempMax-tempMin; end abpMaxTH=mean(tempTH); else [tempMin,jk,tempMax,jl]=seeklocales(abpsig,closeWin,sigLen); abpMaxTH=tempMax-tempMin; end clear j*; clear t*; abpMaxLT=0.4*abpMaxTH; %Seek pulse beats by MinMax method % diffIndex=1; diffIndex=closeWin; %Avoid filter distortion while diffIndex<sigLen tempMin=abpsig(diffIndex); %Initialization tempMax=abpsig(diffIndex); tempIndex=diffIndex; tpeakp=diffIndex; %Avoid initial error tonsetp=diffIndex; %Avoid initial error while tempIndex<sigLen %If no pulses within 2s, then adjust threshold and retry if (tempIndex-diffIndex)>stepWin % tempIndex=diffIndex-closeWin; tempIndex=diffIndex; abpMaxTH=0.6*abpMaxTH; if abpMaxTH<=abpMaxLT abpMaxTH=2.5*abpMaxLT; end break; end if (diff1(tempIndex-1)*diff1(tempIndex+1))<=0 %Candidate fiducial points if (tempIndex+5)<=sigLen jk=tempIndex+5; else jk=sigLen; end if (tempIndex-5)>=1 jj=tempIndex-5; else jj=1; end %Artifacts of oversaturated or signal loss? if (jk-tempIndex)>=5 for ttk=tempIndex:jk if diff1(ttk)~=0 break; end end if ttk==jk break; %Confirm artifacts end end if diff1(jj)<0 %Candidate onset if diff1(jk)>0 [tempMini,tmin,ta,tb]=seeklocales(abpsig,jj,jk); if abs(tmin-tempIndex)<=2 tempMin=tempMini; tonsetp=tmin; end end elseif diff1(jj)>0 %Candidate peak if diff1(jk)<0 [tc,td,tempMaxi,tmax]=seeklocales(abpsig,jj,jk); if abs(tmax-tempIndex)<=2 tempMax=tempMaxi; tpeakp=tmax; end end end if ((tempMax-tempMin)>0.4*abpMaxTH) %evaluation if ((tempMax-tempMin)<2*abpMaxTH) if tpeakp>tonsetp %If more zero-crossing points, further refine! ttempMin=abpsig(tonsetp); ttonsetp=tonsetp; for ttk=tpeakp:-1:(tonsetp+1) if abpsig(ttk)<ttempMin ttempMin=abpsig(ttk); ttonsetp=ttk; end end tempMin=ttempMin; tonsetp=ttonsetp; if peakIndex>0 %If pulse period less than eyeclose, then artifact if (tonsetp-peakp(peakIndex))<(3*closeWin) %too many fiducial points, then reset tempIndex=diffIndex; abpMaxTH=2.5*abpMaxLT; break; end %If pulse period bigger than 2s, then artifact if (tpeakp-peakp(peakIndex))>stepWin peakIndex=peakIndex-1; onsetIndex=onsetIndex-1; if dicroIndex>0 dicroIndex=dicroIndex-1; end end if peakIndex>0 %new pulse beat peakIndex=peakIndex+1; peakp(peakIndex)=tpeakp; onsetIndex=onsetIndex+1; onsetp(onsetIndex)=tonsetp; tf=onsetp(peakIndex)-onsetp(peakIndex-1); to=floor(abpfreq./20); %50ms tff=floor(0.1*tf); if tff<to to=tff; end to=peakp(peakIndex-1)+to; te=floor(abpfreq./2); %500ms tff=floor(0.5*tf); if tff<te te=tff; end te=peakp(peakIndex-1)+te; tff=seekdicrotic(diff1(to:te)); if tff==0 tff=te-peakp(peakIndex-1); tff=floor(tff/3); end dicroIndex=dicroIndex+1; dicron(dicroIndex)=to+tff; tempIndex=tempIndex+closeWin; break; end end if peakIndex==0 %new pulse beat peakIndex=peakIndex+1; peakp(peakIndex)=tpeakp; onsetIndex=onsetIndex+1; onsetp(onsetIndex)=tonsetp; tempIndex=tempIndex+closeWin; break; end end end end end tempIndex=tempIndex+1; %step forwar