Commit 5c4269a9 authored by Kristian Soltesz's avatar Kristian Soltesz
Browse files

Added main code.

parents
function [sys,dia,pp] = dpp(t,y)
% function for compting DeltaPP for arterial pulse pressure signal y(t).
% y (mmHg), t (s). Also returns respiratory angu lar frequency estimate w
%
% sys, dia and pp have the following fields
% k - linear trend
% m - affine offset
% a - amplitude of fit
% b - phase of fit
% w - angular frequency of fit
%
% sys and dia in addition have the field
% J - residual rms per signal rms
% start at t=0
t0 = t-t(1);
% window duration
T = t0(end);
% peak detection
y0 = detrend(y); % linear detrend
mpp = (max(y0)-min(y0))/3; % mean peak prominence
[~,sys.idx] = findpeaks(y0,'MinPeakProminence',mpp);
[~,dia.idx] = findpeaks(-y0,'MinPeakProminence',mpp);
% systolic peaks
sys.y = y(sys.idx);
sys.t = t0(sys.idx);
% diastolic peaks
dia.y = y(dia.idx);
dia.t = t0(dia.idx);
% frequency range
hr = median(1./diff(sys.t)); % median heart rate (beats per second)
wMax = pi*hr/1.5;
wMin = pi/T;
% find frequency
n = 64;
[sys.w,sys.P,sys.wP] = findFrequency(sys.t,sys.y,wMin,wMax,n);
[dia.w,dia.P,dia.wP] = findFrequency(dia.t,dia.y,wMin,wMax,n);
pp.w = mean([sys.w dia.w]);
% fit linear trend, amplitude, and phase
[~,sys.a,sys.b,sys.k,sys.m,sys.J] = fitFrequency(sys.t,sys.y,pp.w);
[~,dia.a,dia.b,dia.k,dia.m,dia.J] = fitFrequency(dia.t,dia.y,pp.w);
% pulse pressure variation
x = sys.a*cos(sys.b)-dia.a*cos(dia.b);
y = sys.a*sin(sys.b)-dia.a*sin(dia.b);
pp.a = sqrt(x^2+y^2);
pp.b = atan2(y,x);
pp.k = sys.k-dia.k;
pp.m = sys.m-dia.m;
% mean pulse pressure
pp.mean = pp.a/(T*pp.w)*(cos(pp.b)-cos(pp.b+T*pp.w))+pp.k*T/2+pp.m;
% minimum and maximal pulse pressure
t0 = mod((-pp.b+acos(-pp.k/(pp.a*pp.w))*[-1 1])/pp.w,2*pi/pp.w); % smallest two t>=0 at which local extrema
tT = t0 + floor((T-t0)/(2*pi/pp.w))*2*pi/pp.w; % largest two t<=T at which local extrema
t = unique([0 t0 tT T]); % candidates
t = t(find(t>=0 & t<=T)); % remove candidates outside range
pp.y = pp.a*sin(pp.w*t+pp.b)+pp.k*t+pp.m;
pp.min = min(pp.y);
pp.max = max(pp.y);
% final output
pp.delta = 100*(pp.max-pp.min)/pp.mean;
return
%function [fmax,P,f] = findFrequency(t,y,fMin,fMax,n)
function [w,P,wP] = findFrequency(t,y,wMin,wMax,n)
% finds most prominent frequency f (Hz) in y (mmHg).
% Chooses between n frequencies, evenly spaced between fMin and fMax
% (both included). Samples are given by t (s( and do not have to be evenly
% spaced in time.
p = polyfit(t,y,1);
y = y-polyval(p,t);
wP = [linspace(wMin,wMax,n)];
P = plomb(y,t,wP/(2*pi)); % 2*pi*f = w
[~,idx] = max(P);
w = wP(idx);
return
function [w,a,b,k,m,J] = fitFrequency(t,y,w)
% Computes linear detrending coefficients k,m and subsequently fits
% harmonic to detrended version of signal y. Harmonic is of angular
%frequency w (rad/s). Returns w,phi,k,m defined through the signal
%k*(a*sin(w*t+b))+m. Also returns residual to signal energy ratio
% compute linear trend
p = polyfit(t,y,1);
k = p(1);
m = p(2);
y0 = y-polyval(p,t);
% compute amplitude and phase
x = [cos(w*t) sin(w*t)]\y0;
a = norm(x);
b = atan2(x(1),x(2));
% residual computation
yHat = a*sin(w*t+b);
e = y0-yHat;
J = rms(e);
return
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment