function NightTimeHRVmodel(n,nos,ifplots)

% % Night-time Heart Rate Variability Model including effect of sleep architecture
% % Version: 1.0
% % Copyright (C) 2016 Mateusz Solinski
% % 
% % 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.
% % 
% % Author: Mateusz Solinski
% % Warsaw University of Technology
% % Faculty of Physics
% % solinski@if.pw.edu.pl
% % http://solinski.fizyka.pw.edu.pl/
% % 
% % Method was first proposed in:
% % M. Solinski, J. Gieraltowski, J. Zebrowski
% % Modeling heart rate variability including the effect of sleep stages, 
% % Chaos 26, 023101 (2016).
% % http://dx.doi.org/10.1063/1.4940762
% % 
% % Please cite the above publication when referencing this material, 
% % and also include the standard citation for PhysioNet:
% % Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. 
% % PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals. 
% % Circulation 101(23):e215-e220 
% % [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215]; 2000 (June 13).

% % USAGE
% % Usage example:
% %   NightTimeHRVmodel(n,nos,ifplot)
% % where: n is number of samples of signal (default 25000), nos is number of signals to generate (default 1) and ifplot parameter
% % enable/disable plotting of generated signals (signals are plotted when ifplot=1). 
% % The output as .txt files included RR time seriee and (optionaly) .png plots are storage at a folder “. /RR” placed in the same
% % directory as Matlab script.

% % This script was tested in Matlab software, version R2015b (8.6.0.267246) and version R2013a (8.1.0.604) (date of testing: 2.03.2016)

  if nargin < 1
    n =  25000;
  end

  if nargin < 2
     nos =  1;
  end

  if nargin < 3
    ifplots =  1;
  end

if (n<500 || n>250000)
   error('The length of signal must be greater than 500 and smaller than 250000 samples.')
end   

    for o=1:nos

        %generate hynogram
        [HH,HHchange,HHdur]=getHypnogram(n); 

        %generate stochastic part of model
        [RR, RRzmiana,K]=StochasticPartToRR(HH,HHchange,HHdur,n);

        %generate gaussian part of model (breathing activity)
        RRGauss=GaussPartToRR(length(RR));
        RR=RR+10*RRGauss;
        RRzmiana=cumsum(RRzmiana); 

        % remove discontinuousness
        RR=discontinuousnessRemove(RRzmiana, RR, length(RR));

        % added characteristic shapes for the lowest values of RR
        RR=excerciseEpisodes(HHchange,RRzmiana, RR,HH, length(RR));

        RR=RR(1:n)*1000;

        if(ifplots==1)
            clf;
            RRplot=figure(1)

            subplot(2,1,1)
            plot(RR,'k-'); 
            xlabel('n')
            ylabel('RR [ms]')
            axis([0 inf 250 1800])
            subplot(2,1,2)
            plot(HH,'g-','LineWidth', 1);
            ax = gca;
            ax.YTickLabel  = {'Exc', 'D', '', 'L','' ,'REM', 'W',''};
            axis([0 inf 0 7])
            xlabel('n')
            ylabel('sleep stage')

            if ~isequal(exist('./RR' , 'dir'),7)
                mkdir('./RR/');
            end

            dlmwrite(['./RR/RR' num2str(o) '.txt'],RR)
            print(RRplot,['./RR/RR' num2str(o)],'-dpng')
        end

    end


end


function [HH, HHchange, HHdur]=getHypnogram(n)

 randLenghtofExcPhase = 30;
 HH=3; 
 
 while(1)
  
  %matrixes of probabilty of transitions between each pair of sleep stage
  if(n < 7200)
  LDp = 0.502; LRp = 0.048; LWp = 0.054;   LEp = 0.395; 
  DLp = 0.749; DRp = 0.009; DWp = 0.018;   DEp = 0.224; 
  RLp = 0.761; RDp = 0.000; RWp = 0.116;   REp = 0.123;
  WLp = 0.041; WDp = 0.000; WRp = 0.917;   WEp = 0.041; 
  ELp = 0.467; EDp = 0.039; ERp = 0.313;   EWp = 0.18;
  else
      if(length(HH)<ceil(n/2))
    LDp = 0.502; LRp = 0.048; LWp = 0.054;   LEp = 0.395; 
    DLp = 0.749; DRp = 0.009; DWp = 0.018;   DEp = 0.224; 
    RLp = 0.761; RDp = 0.000; RWp = 0.116;   REp = 0.123;
    WLp = 0.041; WDp = 0.000; WRp = 0.917;   WEp = 0.041; 
    ELp = 0.467; EDp = 0.039; ERp = 0.313;   EWp = 0.18; 
      else
    LDp = 0.268; LRp = 0.081; LWp = 0.093;   LEp = 0.559; 
    DLp = 0.771; DRp = 0.018; DWp = 0.030;   DEp = 0.181; 
    RLp = 0.703; RDp = 0.003; RWp = 0.172;   REp = 0.122;
    WLp = 0.086; WDp = 0.000; WRp = 0.800;   WEp = 0.114; 
    ELp = 0.428; EDp = 0.009; ERp = 0.336;   EWp = 0.227; 
      end
  end
   r=rand();
   
   % draw the sleep stage with random length (from distributions
   % characterized each stage)
   % numbers of sleep stages are as follows:
   % 1 - deep sleep
   % 3 - light sleep
   % 5 - REM stage
   % 6 - wake stage
   % 0 - excercise episode
   switch HH(length(HH))
       case 1
           if(r<DLp)
              HH=[HH ; 3*ones(ceil(20*random('Gamma',0.6290078212294944,29.092799446920807,1,1)),1)];   
           end
           if(r>DLp && r<=DLp+DRp)
              HH=[HH ; 5*ones(ceil(20*random('Gamma',0.612456,4.8293,1,1)),1)];   
           end
           if(r>DLp+DRp && r<=DLp+DRp+DWp)
              HH=[HH ; 6*ones(ceil(20*random('Exponential',0.204816,1,1)),1)];   
           end
           if(r>DLp+DRp+DWp && r<=1)
              HH=[HH ; 0*ones(randLenghtofExcPhase,1)];   
           end           
       case 3
           if(r<LDp)
              HH=[HH ; 1*ones(ceil(20*random('Gamma',0.3177, 39.9838,1,1)),1)];   
           end
           if(r>LDp && r<=LDp + LRp)
              HH=[HH ; 5*ones(ceil(20*random('Gamma',0.612456,4.8293,1,1)),1)];   
           end
           if(r>LDp + LRp && r<=LDp + LRp + LWp)
              HH=[HH ; 6*ones(ceil(20*random('Exponential',0.204816,1,1)),1)];   
           end
           if(r>LDp + LRp + LWp && r<=1)
              HH=[HH ; 0*ones(randLenghtofExcPhase,1)];   
           end  
       case 5
           if(r<RLp)
              HH=[HH ; 1*ones(ceil(20*random('Gamma',0.3177, 39.9838,1,1)),1)];   
           end
           if(r>RLp && r <= RLp + RDp)
              HH=[HH ; 3*ones(ceil(20*random('Gamma',0.6290078212294944,29.092799446920807,1,1)),1)];  
           end
           if(r > RLp + RDp && r <= RLp + RDp + RWp)
              HH=[HH ; 6*ones(ceil(20*random('Exponential',0.204816,1,1)),1)];   
           end
           if(r>RLp + RDp + RWp && r<=1)
              HH=[HH ; 0*ones(randLenghtofExcPhase,1)];   
           end 
       case 6
           if(r<WLp)
              HH=[HH ; 1*ones(ceil(20*random('Gamma',0.3177, 39.9838,1,1)),1)];   
           end
           if(r>WLp && r <= WLp + WDp)
              HH=[HH ; 3*ones(ceil(20*random('Gamma',0.6290078212294944,29.092799446920807,1,1)),1)];  
           end
           if(r > WLp + WDp && r <= WLp + WDp + WRp)
               HH=[HH ; 5*ones(ceil(20*random('Gamma',0.612456,4.8293,1,1)),1)];   
           end
           if(r>WLp + WDp + WRp  && r<=1)
              HH=[HH ; 0*ones(randLenghtofExcPhase,1)];   
           end
       case 0
           if(r<ELp)
              HH=[HH ; 1*ones(ceil(20*random('Gamma',0.3177, 39.9838,1,1)),1)];   
           end
           if(r>ELp && r <= ELp + EDp)
              HH=[HH ; 3*ones(ceil(20*random('Gamma',0.6290078212294944,29.092799446920807,1,1)),1)];  
           end
           if(r >ELp + EDp && r <= ELp + EDp + ERp)
               HH=[HH ; 5*ones(ceil(20*random('Gamma',0.612456,4.8293,1,1)),1)];   
           end
           if(r>ELp + EDp + ERp  && r<=1)
              HH=[HH ; 6*ones(ceil(20*random('Exponential',0.204816,1,1)),1)];   %mozliwe ze bedzie trzeba dodac warunek na zbyt dlugie wartosci
           end
   end
        
       if(length(HH)>n)
           HH=HH(1:n);
           break
       end
   
 end

   % indexes of sleep stage changes
   HHchange=[];
   for i=1:length(HH)-1
    if(HH(i)~=HH(i+1)) HHchange=[HHchange ; i]; end
   end
   
   HHdur=diff([0 ; HHchange]);
   HHdur=[HHdur ; length(HH)-HHchange(end)];
   HHchange=[HHchange ; length(HH)];
end

function [RR, RRzmiana,k]=StochasticPartToRR(HH,HHchange,HHdur,n)

    %initial parameters
    gamma=0;
    bb=0;
    mi=0;
    k=[];
    y=[];
    i=1;
    RR=[];
    RRzmiana=[];
    window=50;
     
    for s=1:length(HHdur)
        
        RROnlyThatStage=[];
        HHCurr=HH(HHchange(s));

 
       % values of gamma, b and mi parameters change after transition to
       % the next sleep stage
       switch HHCurr
           case 0
               gamma=1.8;
               bb=0.75;
               mi=0.75;
           case 1
               gamma=5;
               bb=0.3;
               mi=0.95;
           case 3
               gamma=2.7;
               bb=0.45;
               mi=0.95;
           case 5
               gamma=2.4;
               bb=0.55;
               mi=0.9;
           case 6
               gamma=2.1;
               bb=0.75;
               mi=0.85;
       end
           
       % randomisation of parameters 
       gamma=normrnd(gamma,0.05*gamma);
       bb=normrnd(bb,0.05*bb);
       mi=normrnd(mi,0.05*mi);
       
       while(1)
          
            %i.i.d. random integers 
            k=[k ; k_gen(gamma, n, 1)]; 
         
        if(i-k(i)>0) 
            average=mean(y(i-k(i):i-1).^2);
            mu=0;
            sigma=1;
            u = rand()-0.5;
            b = sigma / sqrt(2);
            laplace = mu - b * sign(u).* log(1- 2* abs(u));
            y=[y ; 0.5*laplace*sqrt(1+bb*average)]; 
        else
            y=[y ; 0];
        end
          
        if(i<window+1) 
            theta=zeros(i,1);
            j=1:i;
            theta=y(j);
            thIDX=find(k(j)+j'-i<0);
        else
            theta=zeros(window,1);
            j=i-window:i;
            theta=y(j);
            thIDX=find(k(j)+j'-i<0);
        end

        theta(thIDX)=0;
        stoch=sum(theta);
       RR=[RR ; mi+0.025*stoch];
       
        RROnlyThatStage=[RROnlyThatStage ; RR(i)];
        i=i+1;
        RRTotalOnlyThatStage=sum(RROnlyThatStage);
%         HHdur(s)
       
        if(RRTotalOnlyThatStage > HHdur(s))
            RRzmiana=[RRzmiana ; length(RROnlyThatStage)];
            %end of current sleep stage
            break;
        end        
       end      
    end   
end



function k = k_gen(gamma, N, numberOfRandoms)

    % x(end)=N/4 (against N) for faster calculation (there are not significant difference) 
    x = 6:1: floor(N/4); 
    xMin = min(x);
    inversePowerLawPDF = ((gamma-1) / xMin) * (x/xMin) .^ (-gamma);
    % get the CDF numerically
    inversePowerLawCDF = cumsum(inversePowerLawPDF);
    % normalize
    inversePowerLawCDF = inversePowerLawCDF / inversePowerLawCDF(end);
	nearestIndex = find(inversePowerLawCDF >= rand);
	k = round(x(nearestIndex(1)));

end

function  RRGauss=GaussPartToRR(N)

    frequency1 = 0.1; %Hz - Mayer
    frequency2 = 0.25; %Hz - Vagus
    sig1= 0.15;
    sig=2;
    sig2=sig1*sig;

    fs=1;

    nFR=N;
    df=1/nFR;

    w1=2*pi*frequency1;
    w2=2*pi*frequency2;
    sigma1=0.01;
    sigma2=0.02;
    c1=2*pi*sigma1;
    c2=2*pi*sigma2;


    w=[];
    iw=1:nFR;
    w=(iw'-1)*2*pi*df;

    dw1=w-w1;
    dw2=w-w2;
    
    % frequency spectrum modelling breathing influence to HRV
    Hw = sig1*exp(-dw1.^2./(2*c1.^2))+ sig2*exp(-dw2.^2./(2*c2.^2));
    Sw= sqrt([(fs/2)*Hw(1:ceil(nFR/2)) ; (fs/2)*Hw(ceil(nFR/2)+1:-1:1)]); 
    
    % phase randomization 
    ph0=2*pi*rand(ceil(nFR/2)-1,1);

    ph=zeros(nFR,1);
    phI=1:ceil(nFR/2)-1;

    ph(phI'+1)=ph0(phI');
    ph(nFR-phI'+1)=-ph0(phI');
    
    SwC=zeros(2*nFR,1);

    i=1:nFR;
    i=i';
    SwC(2*i-1)=Sw(i).*cos(ph(i));
    SwC(2*i)=Sw(i).*sin(ph(i));
    
    %inverse fourier transform from frequency spectrum
    rr=real(ifft(SwC));

    RRGauss=rr(1:N);

end

function RR=discontinuousnessRemove(RRzmiana, RR,n)

    for j=1:length(RRzmiana)
       
        if(RRzmiana(j)>10 && RRzmiana(j)<n-10)
            fbefore=RR(RRzmiana(j)-9:RRzmiana(j));
            fafter=RR(RRzmiana(j)+1:RRzmiana(j)+10);
            fjoin=[fbefore ; fafter];
            
            %moving mean
            idx=0:10;
            idxALL=1:length(fbefore);
            windows=bsxfun(@plus,idx,idxALL');
            
            fmean=mean(fjoin(windows));
            ifmean=1:length(fmean);
            
            ip1=1:length(fbefore);
            p1=polyfit(ip1',fbefore,1);
            line1=p1(2)+p1(1)*ip1;
            ip2=1:length(fafter);
            p2=polyfit(ip2',fafter,1);
            line2=p2(2)+p2(1)*ip2;            
           
            fmeaninterpolation = interp1(ifmean,fmean,1:0.52:length(fmean));
            diff1=fmeaninterpolation(1:length(fmeaninterpolation)/2)-line1;
            diff2=fmeaninterpolation(length(fmeaninterpolation)/2+1:length(fmeaninterpolation))-line2;
            
            RR(RRzmiana(j)-9:RRzmiana(j)) = RR(RRzmiana(j)-9:RRzmiana(j))+ diff1';
            RR(RRzmiana(j)+1:RRzmiana(j)+10) = RR(RRzmiana(j)+1:RRzmiana(j)+10)+ diff2';

        end
        
    end
end

function RR=excerciseEpisodes(HHchange,RRzmiana, RR, HH, n)

    excIndexesP=[];
    excIndexesL=[];

    for i=1:length(HHchange)

        % choose area around excercise episodes (HH==0)
        if(HH(HHchange(i))==0  )

            if(HHchange(i)>30)
                excIndexesP=[excIndexesP ; RRzmiana(i)];
                excIndexesL=[excIndexesL ; RRzmiana(i)-30];
            else
                excIndexesP=[excIndexesP ; RRzmiana(i)];
                excIndexesL=[excIndexesL ; RRzmiana(i)];
            end
        end
    end

    for g=1:length(excIndexesL)

        lista=RR(excIndexesL(g):excIndexesP(g));
        dodzielnika=floor(normrnd(0.70,0.05)*length(lista));
        dzielnik=excIndexesL(g)+dodzielnika;
        y=[lista(1) ; min(lista) ; lista(end)];
        
        % model of left slope
        dataleftX=[excIndexesL(g) ; excIndexesL(g) + floor((dzielnik-excIndexesL(g))/2) ; dzielnik];
        dataleftY=[y(1) ; y(1)-1.6*(y(1)-y(2))/2 ; y(2)];
        interLeft=interp1(dataleftX,dataleftY,excIndexesL(g):dzielnik,'spline');

        % model of right slope
        dataleftX=[dzielnik ; floor(dzielnik + (excIndexesP(g) -dzielnik)/2) ; excIndexesP(g)];
        dataleftY=[y(2) ; y(2)+0.4*(y(3)-y(2))/2 ; y(3)];
        interRight=interp1(dataleftX,dataleftY,dzielnik:excIndexesP(g),'spline');

        % randomization of slopes
        iterAll=[normrnd(1,0.01,length(interLeft),1).*interLeft' ; ...
            normrnd(1,0.02,length(interRight),1).*interRight'];


        r1=floor(normrnd(0.6,0.025)*dodzielnika);
        r2=floor(normrnd(0.75,0.025)*dodzielnika);

        % adding random little step to left slope
        for i=1:dodzielnika
            if(i>= r1 && i<= r2)
                iterAll(i)=1.02*iterAll(i);
            end
        end

        RR(excIndexesL(g):excIndexesL(g)+length(iterAll)-1)=iterAll;

    end

end