close all
clear 
clc

%% Paroxysmal atrial fibrillation (PAF) generator
%
% Copyright (C) 2017  Andrius Petrenas
% Biomedical Engineering Institute, Kaunas University of Technology
%
% Available under the GNU General Public License version 3
% - please see the accompanying file named "LICENSE"
%
% Related publication: Petrenas, A., Marozas, V., Solosenko, A., Kubilius, R., 
% Skibarkiene, J., Oster, J., & Sornmo, L. (2017). Electrocardiogram modeling 
% during paroxysmal atrial fibrillation: Application to the detection of brief 
% episodes. Physiological Measurement, 38(11), 2058–2080. 
% http://doi.org/10.1088/1361-6579/aa9153
% For more information see help simPAF_ECG_generator!!!

% ---Initial parameters---
rrLength = 50;     % A desired ECG signal length (the number of RR intervals) 
APBrate = 0.10;    % Rate of atrial premature beats (APB). A number between 0 and 0.5
onlyRR = 0;        % 1 - only RR intervals are generated, 0 - multilead ECG is generated

medEpis = 15;       % Median duration of an atrial fibrillation (AF) episode
stayInAF = 1-log(2)/medEpis;   % Probability to stay in AF state
AFburden = 0.8;     % AF burden. 0 - the entire signal is sinus rhythm (SR), 1 - the entire signal is AF

noiseType = 0;      % Type of noise. A number from 0 to 4. 0 - no noise added (noise RMS = 0 mV), 
                    % 1 - motion artefacts, 2 - electrode movement artefacts, 3 - baseline wander, 
                    % 4 - mixture of type 1, type 2 and type 3 noises
noiseRMS = 0.02;    % Noise level in milivolts 

realRRon = 1;       % 1 - real RR series are used, 0 - synthetic
realVAon = 1;       % 1 - real ventricular activity is used, 0 - synthetic
realAAon = 1;       % 1 - real atrial activity is used, 0 - synthetic
% Note: cannot select real atrial activity and synthetic ventricular activity

% ---PAF generator--- 
[simPAFdata, initialParameters] = simPAF_ECG_generator(rrLength, realRRon, realVAon, realAAon, AFburden, stayInAF, APBrate, noiseType, noiseRMS, onlyRR);

% ---Returned data and initial parameters---
% Initial parameters
fibFreqz = initialParameters.fibFreqz;      % Dominant fibrillatory frequency
% Data
rr = simPAFdata.rr;                         % RR intervals in miliseconds
multileadECG = simPAFdata.multileadECG;     % 15-lead ECG
multileadVA = simPAFdata.multileadVA;       % 15-lead ventricular activity
multileadAA = simPAFdata.multileadAA;       % 15-lead atrial activity
multileadNoise = simPAFdata.multileadNoise; % 15-lead physiological noise
QRSindex = simPAFdata.QRSindex;             % QRS index. Shows sample number when R peak occurs
targets_SR_AF = simPAFdata.targets_SR_AF;   % Marks SR and AF beats (1 - AF, 0 - SR)
targets_APB = simPAFdata.targets_APB;       % Marks atrial premature beats (1 - APB, 0 - not APB)
pafBoundaries = simPAFdata.pafBoundaries;   % Start and the end of each PAF episode
pafEpisLength = simPAFdata.pafEpisLength;   % Length (in beats) of each PAF episode
ecgLength = simPAFdata.ecgLength;           % ECG length in samples

%% PPG generator
%
% Copyright (C) 2019  Andrius Solosenko, All Rights Reserved
% Biomedical Engineering Institute, Kaunas University of Technology
%
% PPG generator is a freeware. Free for private, non-commercial use and sharing. 
% Reverse engineering is restricted. You can also download functions 
% [] = gen_PPG() and [] = gen_PPGpulse() at 
% https://biomedicine.ktu.edu/#software 
%
% Related publications: Solosenko, A., Petrenas, A., Marozas, V., & Sornmo, L. (2017). 
% Modeling of the photoplethysmogram during atrial fibrillation. Computers in Biology 
% and Medicine, 81, 130–138. http://doi.org/10.1016/j.compbiomed.2016.12.016
% Paliakaite, B., Petrenas, A., Solosenko, A., & Marozas, V. (2020). Photoplethysmogram
% modeling of extreme bradycardia and ventricular tachycardia. In Mediterranean Conference 
% on Medical and Biological Engineering and Computing – MEDICON 2019, IFMBE Proceedings, 76, 
% 1–10. http://doi.org/10.1007/978-3-030-31635-8_141

RR = simPAFdata.rr; % RR intervals in samples
% Note: simPAFdata.rr is a vector of RR intervals generated by PAF ECG generator;
% however, it is possible to use any available RR intervals (real or simulated). 
% PPG generator is tested for simulating sinus rhythm, atrial fibrillation, premature beats, bradycardia and tachycardia episodes.
pulseType = 1;  % 1 - Type 1, 2 - Type 2, 3 - Type 3, 4 - Type 3 bis, 5 - Type 4
Fd = 500; % PPG sampling frequency

% ---PPG generator---
[PPGmodel, PPGpeakIdx, PPGpeakVal] = gen_PPG(RR, pulseType, Fd); 
% PPGmodel - generated PPG signal
% PPGpeakIdx - sample number when pulse peak occurs 
% PPGpeakVal - PPG pulse peak values

%% PPG artifact generator
%
% Copyright (C) 2020  Birute Paliakaite
% Biomedical Engineering Institute, Kaunas University of Technology
%
% Available under the GNU General Public License version 3
% - please see the accompanying file named "LICENSE"
%
% Related publication: Paliakaite, B., Petrenas, A., Solosenko, A., & Marozas, V. (2021). 
% Modeling of artifacts in the wrist photoplethysmogram: Application to the 
% detection of life-threatening arrhythmias. Biomedical Signal Processing and Control, 66, 102421.
% https://doi.org/10.1016/j.bspc.2021.102421

% ---Initial parameters---
load('artifact_param.mat')	% load predefined artifact parameters:
% RMS_shape - shape parameter for gamma distribution of normalized RMS artifact amplitudes, [1x4] vector
% RMS_scale - scale parameter for gamma distribution of normalized RMS artifact amplitudes, [1x4] vector
% slope_m - mean of estimated PSD slopes, [1x4] vector 
% slope_sd - standard deviation of estimated PSD slopes, [1x4] vector 

typ_artifact = [0 0 1 1];	% Artifact types to generate: 1/0 for Generate/Omit [*device displacement* *forearm motion* *hand motion* *poor contact*]
prob = typ_artifact*(1/sum(typ_artifact));  % Probabilities of artifacts 

dur_mu0 = 10;	% Mean duration of artifact-free intervals in seconds (1/lambda0) 
dur_mu = 10;	% Mean duration of artifacts in seconds (1/lambda1 = 1/lambda2 = 1/lambda3 = 1/lambda4 = 1/lambda) 
duration = length(PPGmodel);	% Duration of signal to generate in samples

for n=1:4
    % generate slope of power spectral density corresponding to particular artifact type 
	S(n) = slope_sd(n)*randn(1,1) + slope_m(n); 
end

% ---Generate artifacts--- 
[artifact] = gen_PPG_artifacts(duration, prob, [dur_mu0 [1 1 1 1] * dur_mu], RMS_shape, RMS_scale, S, Fd);

%% Addition of PPG and artifacts 
%
% make pulsatile component zero-mean with unitary RMS to obtain the correct
% ratio of artifact and pulsatile RMSs
PPGmodel = (PPGmodel-mean(PPGmodel))/rms(PPGmodel-mean(PPGmodel));
PPG = PPGmodel+artifact;

%% Plots
% ---plots for PAF generation---
if onlyRR == 1
    figure
    hold on  
    plot(rr)
    plot(targets_SR_AF, 'r')
    hold off
    legend('RR intervals', 'Targets (SR(0) / AF(1))','Location','SouthEast')
    xlabel('RR interval number')
    ylabel('RR, ms')
    ylim([0 2.2]) 

else
    figure('units','normalized','outerposition',[0 0 1 1])
    h1 = subplot(10,1,1:2);
    plot(QRSindex, rr,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[1 0 0],'MarkerSize',4,'Marker','o','Color',[0 0 0])
    ylabel('rr, ms')
    title('RR intervals')
    set(gca, 'XTick', [], 'XColor', [1 1 1]);
    ylim([200 1600])

    h2 = subplot(10,1,3:10);
    for i = 1:15
    	plot(multileadECG(i,:)'- 3*i, 'k')
        hold on
    end
    stairs(QRSindex, 2*targets_SR_AF-1, 'r')
    ylim([-48 2])
    xlim([0 length(multileadECG(1,:))])
    set(gca, 'YTick', [-45:3:-3 -1 1], 'YTickLabel', {'Z', 'Y', 'X', 'V6', 'V5', 'V4', 'V3', 'V2', 'V1', 'aVF','aVL','aVR','III','II','I', 'non-AF', 'AF'})

    ylabel('ECG lead')
    xlabel('Sample number')
    linkaxes([h1,h2],'x')
    title('Multilead ECG')

end

% ---plots for PPG and PPG artifact generation--- 

figure
ax(1)=subplot(311);plot(PPG, 'k'); title('PPG');
ax(2)=subplot(312);plot(PPGmodel, 'k'); title('Pulsatile component'); hold on; plot(PPGpeakIdx, PPGmodel(PPGpeakIdx), 'ro'); hold off;
ax(3)=subplot(313);plot(artifact, 'k'); title('Artifacts'); xlabel('Sample number')
linkaxes(ax,'x')

%% version 1.3.1 of April 29, 2022