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