function y=edr(varargin)
%
% y = edr(data_type,signal,r_peaks,fs,pqoff, jpoff, gain_ecg, channel, show)
%
% ECG-derived Respiratory (EDR) signal computation from given
% single-lead ECG signal based on the signed area under the QRS complex.
%
% Required Parameters:
%
% data_type
% A 1x1 integer specifying the file data_type
% 0 --> if Matlab file
% 1 --> if record in MIT format
%
% signal
% A Nx1 integer array containing the ECG signal in mV (if data_type=0)
% OR a char string containing record name (if data_type=1)
%
% fs
% A 1x1 integer specifying the sampling frequency in hz (for Matlab variables only)
%
% r_peaks
%
% A Mx1 integer array containing locations of r peaks on signal in s
% OR a char string containing the extension of the annotation file
% with r peaks in samples (e.g. "qrs") (if data_type=1)
%
% optional parameters:
%
% gain_ecg
% A 1x1 integer specifying dig_max/phy_max (default=1)
%
% channel
% A 1x1 integer>1 (default=1) indicating ECG channel (if data_type=1)
%
% pqoff
% A 1x1 integer>0 specifying average distance between PQ junction and
% R peak, in samples
%
% jpoff
% A 1x1 integer>0 specifying average distance between R peak and
% J point, in samples
%
% show
% A 1x1 boolean if true, generates a plot of the estimated
% respiration signal (default = 0).
%
%
% output:
%
% y
% A Mx2 integer matrix containing time in seconds and edr
%
% This code was written by Sara Mariani at the Wyss Institute at Harvard
% based on the open-source PhysioNet code edr.c
% (http://www.physionet.org/physiotools/edr/)
% by George Moody
%
% Author: Sara Mariani, 2014
% Last Modified: November 17, 2014
%
% please report bugs/questions at sara.mariani@wyss.harvard.edu
%
% Example - Extract EDR signal from ECG in PhysioNet's Remote server:
% signal='fantasia/f1o02';
% r_peaks='ecg';
% data_type=1;
% channel=2;
% show=1;
% y=edr(data_type,signal,r_peaks,[],[],[],channel,show);
% wfdb2mat('f1o02')
% [~,signal,Fs,~]=rdmat('f1o02m');
% resp=signal(:,1);
% resp=resp-mean(resp);
% resp=resp*200;
% sec=length(resp)/Fs;
% xax=[.25:.25:sec];
% r=interp1(y(:,1), y(:,2), xax, 'spline');
% figure
% plot(xax,r)
% hold on
% plot([1:length(resp)]/Fs,resp,'r')
% legend('edr','respiratory signal')
% xlabel('time (s)')
%
% see also: ecgpuwave, gqrs
%endOfHelp
%Set default pararameter values
inputs={'data_type','signal','r_peaks','fs','pqoff','jpoff', 'gain_ecg', 'channel' ,'show'};
show=0;
Ninputs=length(inputs);
if nargin>Ninputs
error('Too many input arguments');
end
if nargin<3
error('Not enough input arguments');
end
for n=1:nargin
eval([inputs{n} '=varargin{n};']);
end
for n=nargin+1:Ninputs
eval([inputs{n} '=[];']);
end
% check format and obtain all the features I need
if data_type==0 %matlab
if isempty(gain_ecg)
gain_ecg=1;
end
ECGm=signal*gain_ecg;
if isempty(r_peaks)
error('R peaks locations not provided');
else
tqrs=round(r_peaks*fs); %samples where I have the R peak
end
elseif data_type==1 %wfdb record
if isempty(channel)
channel=1;
end
% read the signal
wfdb2mat(signal);
pp=strfind(signal,'/');
if ~isempty(pp)
signal2=signal(pp(end)+1:length(signal));
else signal2=signal;
end
[~,sig,fs]=rdmat([signal2 'm']);
ECGm=sig(:,channel);
if numel(fs)>1
fs=fs(channel);
end
% read the header
siginfo=wfdbdesc(signal);
siginfo=siginfo(:,channel);
gainstring=siginfo.Gain;
sp=strfind(gainstring,' ');
try
gain_ecg=str2num(gainstring(1:sp-1));
catch
gain_ecg=1;
end
if strfind(gainstring(end-1),'m')
gain_ecg=gain_ecg*1000;
end
ECGm=ECGm*gain_ecg;
% read r_peaks if annotation file
if ischar(r_peaks)
[ann,ty]=rdann(signal,r_peaks);
tqrs=ann(ty=='N');
r_peaks=tqrs/fs;
else
tqrs=round(r_peaks*fs); %samples where I have the R peak
end
else error('format data_type must be 0 or 1');
end
% check if signal is upside-down
if mean(ECGm(tqrs))<mean(ECGm)
ECGm=-ECGm;
end
% EDR COMPUTATION
% 1) filter the signal with a moving window of lpflen=25 ms
lpflen=0.025;
lp=round(lpflen*fs);
w=ones(lp+1,1)./(lp+1);
sample=filter(w,1,ECGm);
% correct for the delay of lp/2
sample(1:round(lp/2))=[];
% correct for the initialization
for i=1:round(lp/2)
sample(i)=mean(ECGm(1:i+round(lp/2)));
end
% 2) find the baseline: moving window again of bflen=1 s
bflen=1;
b=round(bflen*fs);
w2=ones(b+1,1)./(b+1);
baseline=filter(w2,1,sample);
% correct for the delay of b/2
baseline(1:round(b/2))=[];
% correct for the initialization
for i=1:round(b/2)
baseline(i)=mean(sample(1:i+round(b/2)));
end
% 3) find average boundaries of QRS interval
if isempty(jpoff)||isempty(pqoff)
[pqoff, jpoff]=boundaries(sample, baseline, tqrs, fs);
end
% now estimate signed area under QRS complex
sb=sample(1:length(baseline))-baseline;
snar=zeros(size(tqrs));
for i=2:length(tqrs)-1
win=sb(tqrs(i)-pqoff:tqrs(i)+jpoff);
snar(i)=sum(win);
end
if tqrs(end)+jpoff>length(sb)
win=sb(tqrs(end)-pqoff:end);
else
win=sb(tqrs(end)-pqoff:tqrs(end)+jpoff);
end
snar(end)=sum(win);
% now start from signed area and estimate edr
xm=0;
xd=0;
xdmax=0;
xc=0;
x=snar;
r=zeros(size(x));
for i=25:length(x)
d=x(i)-xm;
if xc<500
xc=xc+1;
dn=d/xc;
else
dn=d/xc;
if dn>xdmax
dn=xdmax;
elseif dn<-xdmax
dn=-xdmax;
end
end
xm=xm+dn;
xd=xd+abs(dn)-xd/(xc);
if xd<1
xd=1;
end
xdmax=3*xd/(xc);
r(i)=d/xd;
end
y=r*50;
while (max(y)>127 || min(y)<-128)
y(y<-128)=y(y<-128)+255;
y(y>127)=y(y>127)-255;
end
if(show)
scrsz = get(0,'ScreenSize');
figure('Position',...
[0.05*scrsz(3) 0.05*scrsz(4) 0.8*scrsz(3) 0.89*scrsz(4)],...
'Color',[1 1 1]);
ax(1)=subplot(211);
plot([1:length(sample)]/fs,sample);
hold on;
plot([1:length(baseline)]/fs,baseline,'g');
plot((tqrs-pqoff)/fs,mean(ECGm)*ones(size(tqrs)),'*m');
plot((tqrs+jpoff)/fs,mean(ECGm)*ones(size(tqrs)),'*c');
legend('filtered ecg','baseline','window start','window end');
set(gca,'fontsize',18);
xlabel('time (s)','fontsize',18);
ylim([mean(ECGm)-5*std(ECGm) mean(ECGm)+5*std(ECGm)]);
ax(2)=subplot(212);
plot(r_peaks,y,'r');
title('edr','fontsize',18);
set(gca,'fontsize',18);
xlabel('time (s)','fontsize',18);
linkaxes(ax,'x');
end
y=[r_peaks y];
end
%%%% Helper function %%%%%%
function[pqoff, jpoff]=boundaries(sample, baseline, tqrs, fs)
% estimate the noise level
sb=sample(1:length(baseline))-baseline;
nlest=mean(abs(sb));
display(['The estimated noise level is ' num2str(nlest) ' microvolts']);
dlthresh=2*nlest;
dlthmax=1200;
dlthmin=140;
if dlthresh>dlthmax, dlthresh=dlthmax;
elseif dlthresh<dlthmin, dlthresh=dlthmin;
end
% determine if samples are baseline
vwindow=100;
twin1=0.033;
twin2=0.067;
% time of the 51st QRS
last=tqrs(51);
sample2=sample(1:last);
bline=zeros(size(sample2));
% a sample is baseline if I have twin1 or twin2 consecutive samples
% that vary in amplitude by no more than dlthresh
for i=1:length(sample2)-twin1*fs
vmax=sample(i);
vmin=sample(i);
if abs(baseline(i)-vmax)<vwindow, twindow=twin1;
else twindow=twin2;
end
ww=sample(i:i+round(twindow*fs));
if max(ww)-min(ww)<dlthresh
bline(i)=1;
end
end
% for first 50 beats, look for PQ junction and J point
tlim2=0.060;
tlim3=0.100;
PQ=zeros(50,1);
J=zeros(50,1);
for j=1:50
% search to the left
try
w=bline(round(tqrs(j)-tlim2*fs):tqrs(j)-1);
catch
display(j);
w=bline(1:tqrs(j)-1);
end
f=find(w);
if numel(f)>0
PQ(j)=length(w)-max(f)+1;
else
PQ(j)=length(w);
end
% search to the right
w=bline(tqrs(j)+1:round(tqrs(j)+tlim3*fs));
f=find(w);
if numel(f)>0
J(j)=min(f);
else
J(j)=length(w);
end
end
% incremental average
pqoff=PQ(1);
for i=1:length(PQ)
if PQ(i)<pqoff
pqoff=pqoff-1;
elseif PQ(i)>pqoff
pqoff=pqoff+1;
end
end
jpoff=J(1);
for i=1:length(J)
if J(i)<jpoff
jpoff=jpoff-1;
elseif J(i)>jpoff
jpoff=jpoff+1;
end
end
end