function [QRSindex, rr, multileadVA, ecgLength] = simPAF_gen_multilead_VA(ecgLength, targets_SR_AF, rr, AFburden, realVAon, realAAon)
%
% [] = simPAF_gen_multilead_VA() returns multilead (15 lead) ventricular
% activity. A set of 100 15-lead ECGs with SR selected from the PTB Diagnostic
% ECG Database is used as a basis for modeling ventricular activity. The ECGs
% of the PTB database are first subjected to baseline removal and QRST delineation.
% The original T-waves are then resampled to a fixed width and, depending on
% the type of rhythm, width-adjusted to match prevailing heart rate. Since
% the original ECGs last just for about 2 min, QRST complexes are subjected
% to repeated concatenation until desired length of ECG is obtained. The TQ
% interval is interpolated using a cubic spline interpolation.
%
% 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"
%
% Generated leads:
% multileadVA(1,:) - I multileadVA(7,:) - V1 multileadVA(13,:) - X
% multileadVA(2,:) - II multileadVA(8,:) - V2 multileadVA(14,:) - Y
% multileadVA(3,:) - III multileadVA(9,:) - V3 multileadVA(15,:) - Z
% multileadVA(4,:) - aVR multileadVA(10,:) - V4
% multileadVA(5,:) - aVL multileadVA(11,:) - V5
% multileadVA(6,:) - aVF multileadVA(12,:) - V6
rr = round(rr*1000); % rr intervals in miliseconds
% Generate PQRST complexes
switch realVAon
case 0 % Generate synthetic PQRST complexes
[data.pqrst, data.Qind, data.Rind, data.Sind, data.Tind] = simPAF_gen_syn_VA(100);
case 1 % Load random patch of PQRST complexes extracted from PTB database
sigNum = randi([1 100]);
load('DATA_PQRST_real')
data.pqrst = DATApqrst(sigNum).pqrst;
data.Qind = DATApqrst(sigNum).Qind;
data.Rind = DATApqrst(sigNum).Rind;
data.Sind = DATApqrst(sigNum).Sind;
data.Tind = DATApqrst(sigNum).Tind;
end
% Original PQRST is subjected to repeated concatenation until required
% number of beats is obtained
arraySize = length(data.pqrst(1,:,1));
pqrstLength = length(data.pqrst(1,1,:));
pqrst = zeros(15, ecgLength,pqrstLength);
j = 1;
for i = 1:ecgLength+1 % Additonal is required to meat exact number of RR intervals
pqrst(:,i,:) = data.pqrst(:,j,:);
j = j + 1;
if j > arraySize
j = 1;
end
end
%Load QRST complexes for all leads
cnt = 0;
for lead = 1:15
ecgSig = [];
AFend = 0;
cnt = cnt + 1;
switch realAAon
case 0 % Synthetic atrial activity
rIndex = data.Rind - data.Qind;
Qind = data.Qind; Rind = data.Rind - data.Qind;
Sind = data.Sind - data.Qind; Tind = data.Tind - data.Qind;
for beatNr = 1:ecgLength
QRSTtemp = pqrst(lead,beatNr,Qind+1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
QRSTtempNext = pqrst(lead,beatNr+1,Qind+1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
% Divide PQRST into two parts
QRST_PS = QRSTtemp(1,1:Sind); % The PQRS part
QRST_ST = QRSTtemp(1,Sind+1:Tind); % The ST part
% Resample T wave according to current RR interval
% Signal is prolonged in order to avoid boundary effects of resampling
QRST_ST = [QRST_ST(1,1)*ones(1,20) QRST_ST QRST_ST(1,end)*ones(1,20)];
% Gradual prolongation of T wave after AF is terminated
if targets_SR_AF(beatNr) == 0
if AFend == 1
ST_length = round(length(QRST_ST)*sqrt((0.3+m*0.1)*rr(beatNr)/1000));
m = m + 1;
if m == 7
AFend = 0;
end
else
ST_length = round(length(QRST_ST)*sqrt(rr(beatNr)/1000));
end
else
ST_length = round(length(QRST_ST)*sqrt(0.35));
AFend = 1;
m = 1;
end
QRST_ST = resample(QRST_ST, ST_length+40, length(QRST_ST)); % T wave length correction QT = QTc*sqrt(RR) ~ T = Tc*sqrt(RR)
QRSTc = [QRST_PS QRST_ST(21:end-20)]; % Corrected PQRST
% Find TQ length
TQlength = rr(beatNr) - length(QRSTc);
% Protect against the error of to low heart rate
if TQlength < 2
TQlength = 2;
end
% Interpolate TQ interval
x = [1 TQlength];
xi = 1:1:TQlength;
y = [QRSTc(1,end) QRSTtempNext(1, 1)];
TQ = interp1(x,y,xi,'linear');
ecgSig = [ecgSig QRSTc TQ];
if lead == 15
rIndex = [rIndex (rIndex(1,end) + TQlength + length(QRSTc))];
end
end
case 1 % Real atrial activity
if AFburden == 1 % Entire signal is AF
rIndex = data.Rind - data.Qind;
else
rIndex = data.Rind;
end
for beatNr = 1:ecgLength
if beatNr == 1
if AFburden == 1
Qind = data.Qind; Rind = data.Rind - data.Qind;
Sind = data.Sind - data.Qind; Tind = data.Tind - data.Qind;
QRSTtemp = pqrst(lead,beatNr,Qind + 1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
else
Qind = data.Qind; Rind = data.Rind;
Sind = data.Sind; Tind = data.Tind;
QRSTtemp = pqrst(lead,beatNr,1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
end
else
if targets_SR_AF(beatNr) == 0
if targets_SR_AF(beatNr-1) == 0
QRSTtemp = pqrst(lead,beatNr,1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
Qind = data.Qind; Rind = data.Rind;
Sind = data.Sind; Tind = data.Tind;
else % targets_SR_AF(beatNr-1) == 1
QRSTtemp = pqrst(lead,beatNr,Qind+1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
Qind = data.Qind; Rind = data.Rind - data.Qind;
Sind = data.Sind - data.Qind; Tind = data.Tind - data.Qind;
end
else % targets_SR_AF(beatNr) == 1
if targets_SR_AF(beatNr-1) == 1
QRSTtemp = pqrst(lead,beatNr,Qind+1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
Qind = data.Qind; Rind = data.Rind - data.Qind;
Sind = data.Sind - data.Qind; Tind = data.Tind - data.Qind;
else % targets_SR_AF(beatNr-1) == 0
QRSTtemp = pqrst(lead,beatNr,1:end);
QRSTtemp = simPAF_correct_baseline(QRSTtemp);
Qind = data.Qind; Rind = data.Rind;
Sind = data.Sind; Tind = data.Tind;
end
end
end
if beatNr == length(targets_SR_AF)
if targets_SR_AF(beatNr) == 0
RindNext = data.Rind;
QRSTtempNext = pqrst(lead,beatNr,1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
else
RindNext = data.Rind-data.Qind;
QRSTtempNext = pqrst(lead,beatNr,Qind+1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
end
else
if targets_SR_AF(beatNr+1) == 0
if targets_SR_AF(beatNr) == 1
RindNext = data.Rind-data.Qind;
QRSTtempNext = pqrst(lead,beatNr+1,Qind+1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
else
RindNext = data.Rind;
QRSTtempNext = pqrst(lead,beatNr+1,1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
end
else
if targets_SR_AF(beatNr) == 0
RindNext = data.Rind;
QRSTtempNext = pqrst(lead,beatNr+1,1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
else
RindNext = data.Rind-data.Qind;
QRSTtempNext = pqrst(lead,beatNr+1,Qind+1:end);
QRSTtempNext = simPAF_correct_baseline(QRSTtempNext);
end
end
end
% Divide PQRST into two parts
QRST_PS = QRSTtemp(1,1:Sind); % The PQRS part
QRST_ST = QRSTtemp(1,Sind+1:Tind); % The ST part
% Resample T wave according to current RR interval
% Signal is prolonged in order to avoid boundary effects of resampling
QRST_ST = [QRST_ST(1,1)*ones(1,20) QRST_ST QRST_ST(1,end)*ones(1,20)];
% Gradual prolongation of T wave after AF is terminated
if targets_SR_AF(beatNr) == 0
if AFend == 1
ST_length = round(length(QRST_ST)*sqrt((0.3+m*0.1)*rr(beatNr)/1000));
m = m + 1;
if m == 7
AFend = 0;
end
else
ST_length = round(length(QRST_ST)*sqrt(rr(beatNr)/1000));
end
else
ST_length = round(length(QRST_ST)*sqrt(0.35));
AFend = 1;
m = 1;
end
QRST_ST = resample(QRST_ST, ST_length+40, length(QRST_ST)); % T wave length correction QT = QTc*sqrt(RR) ~ T = Tc*sqrt(RR)
QRSTc = [QRST_PS QRST_ST(21:end-20)]; % Corrected PQRST
% Find TQ length
TQlength = rr(beatNr) - length(QRSTc) + Rind - RindNext;
% Protect against the error of to low heart rate
if TQlength < 2
TQlength = 2;
end
% Interpolate TQ interval
x = [1 TQlength];
xi = 1:1:TQlength;
y = [QRSTc(1,end) QRSTtempNext(1, 1)];
TQ = interp1(x,y,xi,'linear');
ecgSig = [ecgSig QRSTc TQ];
if lead == 15
rIndex = [rIndex (rIndex(1,end) + TQlength + length(QRSTc) - Rind + RindNext)];
end
end
if cnt == 1
multileadVA = zeros(15, length(ecgSig));
end
end
if lead == 15
rr = diff(rIndex);
QRSindex = rIndex(1:end-1);
end
multileadVA(lead, :) = ecgSig;
end
ecgLength = length(ecgSig);
end