function [ N_STATES, AS, VS, PS, QS, z ] = ecggen( N, hrmean )
% Produces synthetic ECG with the following outputs:
% s: ECG (mV)
% ipeaks: labels for PQRST peaks: P(1), Q(2), R(3), S(4), T(5)
% A zero lablel is output otherwise ... use R=find(ipeaks==3); 
% to find the R peaks s(R), etc. 
% 
% Operation uses the following parameters (default values in []s):
% sfecg: ECG sampling frequency [256 Hertz]
% N: approximate number of heart beats [256]
% Anoise: Additive uniformly distributed measurement noise [0 mV]
% hrmean: Mean heart rate [60 beats per minute]
% hrstd: Standard deviation of heart rate [1 beat per minute]
% lfhfratio: LF/HF ratio [0.5]
% sfint: Internal sampling frequency [256 Hertz]
% Order of extrema: [P Q R S T]
% ti = angles of extrema [-70 -15 0 15 100] degrees
% ai = z-position of extrema [1.2 -5 30 -7.5 0.75]
% bi = Gaussian width of peaks [0.25 0.1 0.1 0.1 0.4]
% Copyright (c) 2003 by Patrick McSharry & Gari Clifford, All Rights Reserved  
% See IEEE Transactions On Biomedical Engineering, 50(3), 289-294, March 2003.
% Contact P. McSharry (patrick@mcsharry.net) or G. Clifford (gari@mit.edu)

%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
% 
% ecgsyn.m and its dependents are freely availble from Physionet - 
% http://www.physionet.org/ - please report any bugs to the authors above.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ECG generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set parameter default values
sfecg = 512;
%N = 60;
Anoise = 0;
%hrmean = 60;
hrstd = 1;
lfhfratio = 0.5;
sfint = 512;

%     P  Q  R  S  T  
ti = [-22 -15 0 15 100];

% convert to radians
ti = ti*pi/180;
%  P   Q   R   S   T   
ai = [1.2 -5 30 -7.5 0.75];

% Gaussian width of each attractor
%    P    Q   R   S   T   
bi = [0.15 0.1 0.1 0.1 0.4];

% adjust extrema parameters for mean heart rate 
hrfact =  sqrt(hrmean/60);
hrfact2 = sqrt(hrfact);
bi = hrfact*bi;
ti = [hrfact2 hrfact 1 hrfact hrfact2].*ti;

% Paramters for atrium and ventricle
ERP_A=0.15; % in s
ERP_V=0.18; % in s
RRP_A=0.15; % in s
RRP_V=0.18; % in s

% check that sfint is an integer multiple of sfecg 
q = round(sfint/sfecg);
qd = sfint/sfecg;
if q ~= qd
   error(['Internal sampling frequency (sfint) must be an integer multiple ' ... 
'of the ECG sampling frequency (sfecg). Your current choices are: ' ... 
'sfecg = ' int2str(sfecg) ' and sfint = ' int2str(sfint) '.']);
end

% define frequency parameters for rr process 
% flo and fhi correspond to the Mayer waves and respiratory rate respectively
flo = 0.1;
fhi = 0.25;
flostd = 0.01;
fhistd = 0.01;

fid = 1;
fprintf(fid,'ECG sampled at %d Hz\n',sfecg);
fprintf(fid,'Approximate number of heart beats: %d\n',N);
fprintf(fid,'Measurement noise amplitude: %d \n',Anoise);
fprintf(fid,'Heart rate mean: %d bpm\n',hrmean);
fprintf(fid,'Heart rate std: %d bpm\n',hrstd);
fprintf(fid,'LF/HF ratio: %g\n',lfhfratio);
fprintf(fid,'Internal sampling frequency: %g\n',sfint);
fprintf(fid,'      P  Q  R  S  T\n'); 
fprintf(fid,'ti = [%g %g %g %g %g] radians\n',ti(1),ti(2),ti(3),ti(4),ti(5));
fprintf(fid,'ai = [%g %g %g %g %g]\n',ai(1),ai(2),ai(3),ai(4),ai(5));
fprintf(fid,'bi = [%g %g %g %g %g]\n',bi(1),bi(2),bi(3),bi(4),bi(5));


% calculate time scales for rr and total output
sampfreqrr = 1;
trr = 1/sampfreqrr; 
%tstep = 1/sfecg;
rrmean = (60/hrmean);	 
Nrr = 2^(ceil(log2(N*rrmean/trr)));

% compute rr process
rr0 = rrprocess(flo,fhi,flostd,fhistd,lfhfratio,hrmean,hrstd,sampfreqrr,Nrr);

% upsample rr time series from 1 Hz to sfint Hz
rr = interp(rr0,sfint);

% make the rrn time series
dt = 1/sfint;
rrn = zeros(length(rr),1);
tecg=0;
i = 1;
while i <= length(rr)
   tecg = tecg+rr(i);
   ip = round(tecg/dt);
   rrn(i:ip) = rr(i);
   i = ip+1;
end 
Nt = ip;

% integrate system using fourth order Runge-Kutta
fprintf(fid,'Integrating dynamical system\n');
x0 = [1,0,0.04];
Tspan = [0:dt:(Nt-1)*dt];
[T,X0] = ode45('derivsecgsyn',Tspan,x0,[],rrn,sfint,ti,ai,bi);

% downsample to required sfecg
X = X0(1:q:end,:);

% Total number of states
N_STATES=size(X,1);

fprintf(fid,'N states: %d\n',N_STATES);

% extract R-peaks times
ipeaks = detectpeaks(X, ti, sfecg);

% Scale signal to lie between -0.4 and 1.2 mV
z = X(:,3);
zmin = min(z);
zmax = max(z);
zrange = zmax - zmin;
z = (z - zmin)*(1.6)/zrange -0.4;

% Discretize the signal
P=find(ipeaks==1);
Q=find(ipeaks==2);
R=find(ipeaks==3);
S=find(ipeaks==4);
T=find(ipeaks==5);



% Atrium signal generation
% Beginning of the atrium signal for pacemaker
% Divide in two the width of the gaussian for the P signal and then divide by tstep 
Ab=sfecg*bi(1)/2;

% Make it integer
Ab=floor(Ab);

% Compute the length of the ERP and RRP for atrium
len_ERP_A=ERP_A/dt;
len_ERP_A=floor(len_ERP_A);

len_RRP_A=RRP_A/dt;
len_RRP_A=floor(len_RRP_A);

% Generate the P signal
PS=zeros(size(z));

% Matrix for the indices of the AS signal
AS=zeros(size(P));

for i=1:size(P,1)
    AS(i)=P(i)-Ab;
    for j_erp=1:len_ERP_A,
        if P(i)-Ab+j_erp-1>size(PS,1)
            break;
        end
        PS(P(i)-Ab+j_erp-1)=1;
    end
    for j_rrp=1:len_RRP_A,
        if P(i)-Ab+j_erp-1+j_rrp>size(PS,1)
            break;
        end
        PS(P(i)-Ab+j_erp-1+j_rrp)=2;
    end    
end

% Ventricle signal generation
% Generate the Q signal
QS=zeros(size(z));

% Compute the length of the ERP and RRP for atrium
len_ERP_V=ERP_V/dt;
len_ERP_V=floor(len_ERP_V);

len_RRP_V=RRP_V/dt;
len_RRP_V=floor(len_RRP_V);

% Matrix for the indices of the VS signal
VS=zeros(size(P));

for i=1:size(Q,1)
    VS(i)=Q(i);
    
    for j_erp=1:len_ERP_V,
        if Q(i)+j_erp-1>size(QS,1)
            break;
        end
        QS(Q(i)+j_erp-1)=1;
    end
    for j_rrp=1:len_RRP_V,
        if Q(i)+j_erp-1+j_rrp>size(QS,1)
            break;
        end
        QS(Q(i)+j_erp-1+j_rrp)=2;
    end    

    
    QS(Q(i))=1;
end

%periods=zeros(size(z));

%cnti=0;
%for i=1:size(rr0,1),
%    cnti=cnti+rr0(i);
%end

%cnt=0;
%for i=1:size(rr0,1),
%    tt=floor(rr0(i)*sfecg);
%    cnt=cnt+tt;
    
%    if cnt>size(z,1),
%        break;
%    end
%    periods(cnt)=1;
%end

end

