function [pqrst, Qind, Rind, Sind, Tind] = simPAF_gen_syn_VA(num)
%
% [] = simPAF_gen_syn_VA() returns simulated QRST complexes. 
% Ventricular activity is simulated using a dynamic ECG model proposed in 
% Sameni et al. 2007, which is an extended version of a single channel
% ECG simulator originally proposed in McSharry et al. 2003. By using this
% model, three orthogonal Frank leads of N samples length are generated.
% Ventricular activity is generated in 3 Frank leads and then 12-lead signals
% are obtained by applying a transformation matrix.
%
% 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"
%
N = 1000;                       % Signal length in samples 
fs = 1000;                      % sampling rate
% num - number of PQRST compelxes
pqrst = zeros(15,num,560);
pqrstFrank = zeros(3,560);

% Dipole parameters 
F = 1;                          % heart rate 
R0 = simPAF_Rotate3D(0,0,0);    % dipole rotation matrices (tetax,tetay,tetaz) 
Lambda = eye(3); 
teta0 = -pi/2;                  % initial phase of the ECG 
 
%% QRS
Qw = simPAF_gen_rand_num(0.05, 0.08);
Rw = simPAF_gen_rand_num(0.05, 0.08);
Sw = simPAF_gen_rand_num(0.05, 0.08);

alphaiQRS.x = [ simPAF_gen_rand_num(-0.4, -0.05)  simPAF_gen_rand_num(0.4, 1.5)     0]; 
biQRS.x     = [Qw    Rw    Sw]; 
tetaiQRS.x  = [-0.1    0    0.1]; 

alphaiQRS.y = [ 0    simPAF_gen_rand_num(0.1, 0.7)   simPAF_gen_rand_num(-0.3, -0.05)]; 
biQRS.y     = [Qw   Rw   Sw]; 
tetaiQRS.y  = [-0.1   0   0.1]; 
  
alphaiQRS.z = [simPAF_gen_rand_num(-0.4, -0.05)   0  simPAF_gen_rand_num(0.1, 1)]; 
biQRS.z     = [Qw   Rw  Sw]; 
tetaiQRS.z  = [ -0.1   0   0.1];

%% T wave
Tw = simPAF_gen_rand_num(0.5, 0.7);

Txa = simPAF_gen_rand_num(0.02, 0.08);
Tya = simPAF_gen_rand_num(0.01, 0.03);
Tza = simPAF_gen_rand_num(-0.02, -0.06);

alphaiT.x = [Txa   2*Txa   3*Txa]; 
biT.x     = [Tw     Tw/2    Tw/4]; 
tetaiT.x  = [1.1     1.4    1.6];

alphaiT.y = [Tya   2*Tya   3*Tya]; 
biT.y     = [Tw     Tw/2    Tw/4]; 
tetaiT.y  = [1.1     1.4    1.6];

alphaiT.z = [Tza   2*Tza   3*Tza]; 
biT.z     = [Tw     Tw/2    Tw/4]; 
tetaiT.z  = [1.1     1.2    1.4];
%%
alphai.x = [alphaiQRS.x alphaiT.x];
alphai.y = [alphaiQRS.y alphaiT.y];
alphai.z = [alphaiQRS.z alphaiT.z];

bi.x = [biQRS.x biT.x];
bi.y = [biQRS.y biT.y];
bi.z = [biQRS.z biT.z];

tetai.x = [tetaiQRS.x tetaiT.x];
tetai.y = [tetaiQRS.y tetaiT.y];
tetai.z = [tetaiQRS.z tetaiT.z];

f = randi([10 30])/10;
% Generate phase for each variability related frequency
if rand(1,1) > 0.5
    phaseX = 0;
else 
    phaseX = pi;
end

if rand(1,1) > 0.5
    phaseY = 0;
else 
    phaseY = pi;
end

if rand(1,1) > 0.5
    phaseZ = 0;
else 
    phaseZ = pi;
end

%
for n = 1:num
    for m = 1:6   
        rangeLow = alphai.x(1,m) - alphai.x(1,m)/10;
        rangeHigh = alphai.x(1,m) + alphai.x(1,m)/10;
        alphaiTemp.x(1,m) = simPAF_gen_QRST_morf_variab(alphai.x(1,m), f, phaseX, rangeLow, rangeHigh, n);    
        
        rangeLow = alphai.y(1,m) - alphai.y(1,m)/10;
        rangeHigh = alphai.y(1,m) + alphai.y(1,m)/10;
        alphaiTemp.y(1,m) = simPAF_gen_QRST_morf_variab(alphai.y(1,m), f, phaseY, rangeLow, rangeHigh, n);
        
        rangeLow = alphai.z(1,m) - alphai.z(1,m)/10;
        rangeHigh = alphai.z(1,m) + alphai.z(1,m)/10;       
        alphaiTemp.z(1,m) = simPAF_gen_QRST_morf_variab(alphai.z(1,m), f, phaseZ, rangeLow, rangeHigh, n);
    end
    
       
[DIP, ~] = simPAF_DipoleGenerator(N,fs,F,alphaiTemp,bi,tetai,teta0); 
pqrstFrank = R0*Lambda*[DIP.x(51:610) ; DIP.y(51:610) ; DIP.z(51:610)]; 

% Dower transform
D = [  
-0.515  0.157   -0.917;
 0.044  0.164   -1.387;
 0.882  0.098   -1.277;
 1.213  0.127   -0.601;
 1.125  0.127   -0.086;
 0.831  0.076    0.230;
 0.632 -0.235    0.059;
 0.235  1.066   -0.132];
  
 pqrst_v1_v6_I_II = D*pqrstFrank; 

 pqrst(1:8,n,:) = pqrst_v1_v6_I_II(1:8,:); % V1, V2, V3, V4, V5, V6, I, II, 
 pqrst(9,n,:) = pqrst_v1_v6_I_II(8,:) - pqrst_v1_v6_I_II(7,:);     % III
 pqrst(10,n,:) = -(pqrst_v1_v6_I_II(7,:) + pqrst_v1_v6_I_II(8,:))/2; % aVR
 pqrst(11,n,:) =  pqrst_v1_v6_I_II(7,:) - pqrst_v1_v6_I_II(8,:)/2;   % aVL
 pqrst(12,n,:) =  pqrst_v1_v6_I_II(8,:) - pqrst_v1_v6_I_II(7,:)/2;   % aVF
    
 pqrst(13,n,:) = pqrstFrank(1,:); % X
 pqrst(14,n,:) = pqrstFrank(2,:); % Y
 pqrst(15,n,:) = pqrstFrank(3,:); % Z
end
 
 [~,qrsMaxInd] = max(pqrstFrank(1,1:300));
 
Qind = qrsMaxInd-50;
Rind = qrsMaxInd;
Sind = qrsMaxInd+50;
Tind = length(pqrstFrank(1,:));

pqrst_temp = pqrst;
pqrst(1:6,:,:) = pqrst(7:12,:,:);
pqrst(7:12,:,:) = pqrst_temp(1:6,:,:);
end