%% Extract Rhythm parameters form ECG
%  This method will try a different number of methods to find the QRS
%  complexes in case one fails
%
%  Inputs
%   -   ecg     The ECG signal
%   -   fs      The sampling frequency in Hz
%   -   QRS     The QRS complex locations (optional)
%  Outputs
%   -   parameters      This method computes 44 parameters and returns them
%                       in a matrix
%   -   parameterNames  Cell array with names of the parameters
%
%  Pietro Bonizzi
%  Joel Karel
%  Maastricht University
%  2017
function [parameters,parameterNames] = extractParameters(ecg,fs,QRS)

%% enter parameter names
%{
parameterNames = {'P_LF'; 'f_LF'; 'P_HF'; 'f_HF'; 'meanRR'; 'SDRR'; 'rMSSDRR'; 'pNN10'; 'pNN20'; ...
    'pNN30'; 'pNN40'; 'pNN50'; 'pNN60'; 'pNN70'; 'pNN80';'pNN90'; 'pNN100'; 'saen1'; 'saen2'; ...
    'V0'; 'V1'; 'LV2'; 'UV2'; 'numbercomponents95percQRST'; 'SimQRST1'; 'SimQRST2'; ...
    'SimQRST3'; 'SimQRST4'; 'SimQRST5'; 'numbercomponents95percP'; 'SimP1'; 'SimP2'; ...
    'SimP3'; 'SimP4'; 'SimP5'; 'SCN'; 'DF'; 'fWP'; 'fWPMAW'; 'rhe';  'SampEn'; 'SampEnMAW'; ... 
    'OI'; 'SE'; 'FWA'; 'differenceMediaRRmeanRR'};
%}
%saen1 and saen2 removed
parameterNames = {'P_LF'; 'f_LF'; 'P_HF'; 'f_HF'; 'meanRR'; 'SDRR'; 'rMSSDRR'; 'pNN10'; 'pNN20'; ...
    'pNN30'; 'pNN40'; 'pNN50'; 'pNN60'; 'pNN70'; 'pNN80';'pNN90'; 'pNN100'; ...
    'V0'; 'V1'; 'LV2'; 'UV2'; 'numbercomponents95percQRST'; 'SimQRST1'; 'SimQRST2'; ...
    'SimQRST3'; 'SimQRST4'; 'SimQRST5'; 'numbercomponents95percP'; 'SimP1'; 'SimP2'; ...
    'SimP3'; 'SimP4'; 'SimP5'; 'SCN'; 'DF'; 'fWP'; 'fWPMAW'; 'rhe';  'SampEn'; 'SampEnMAW'; ... 
    'OI'; 'SE'; 'FWA'; 'differenceMediaRRmeanRR'};

%% check whether QRS complexes are provided and otherwise compute them
if(nargin<3)
    QRS = findQRScomplexes(ecg,fs);
end
    
if length(QRS)< length(ecg)/fs*0.4 % 0.66
    parameters=inf*ones(1,44);
    return
end

% check noisy signal
differenceMediaRRmeanRR = abs(median(diff(QRS))-mean(diff(QRS)))/fs;

%% RR-based measures

%% 1) spectral indices
% interpolate to obtain evenly sampled RR series at 4 Hz
RR = diff(QRS)/fs;
y = interp1(QRS(2:end)/fs,RR,QRS(2)/fs:1/4:QRS(end)/fs,'spline');
[B,A] = cheby2(3,20,[0.05/(4/2) 0.6/(4/2)]);
try
y = filtfilt(B,A,y);
catch
y = padarray(y,[0 length(y)],'pre');
y = filtfilt(B,A,y);
end

% find optimal AR model for resampled RR series
maxlag = 100;
AIC = zeros(maxlag,1);
for p = 1:maxlag
    clear arparameters
    if p<length(y)
    arparameters = arburg(y,p);
    arparameters(1) = [];
    else 
        break
    end
    matrixA = zeros(length(y)-p,p);
    for i = p+1:length(y)
        matrixA(i-p,:) = y(i-1:-1:i-p)';
    end
    y_fit = matrixA*arparameters';
    AIC(p) = length(y)*log(1/length(y)*sum((y(p+1:end)-y_fit').^2))+2*(p);
    if p>2
        if AIC(p-2)>AIC(p-1) && AIC(p)>AIC(p-1)
            break
        end
    end
end
if p == 100
    [~, optimalmodelorder] = min(AIC);
else
optimalmodelorder = p-1;
end

condition = 1;
while condition
    % find LF and HF frequencies and powers
    [PSD,F]=pburg(y,optimalmodelorder,[],4,'onesided');
    
    intLF = find(F>=0.03 & F<=0.15);
    [peaks,idxpeaks] = findpeaks(PSD);
    tempidxLF = find(ismember(idxpeaks,intLF));
    [P_LF,idxLF] = max(peaks(tempidxLF));
    f_LF = F(idxpeaks(tempidxLF(idxLF)));
    
    intHF = find(F>0.15 & F<=0.40);
    tempidxHF = find(ismember(idxpeaks,intHF));
    [P_HF,idxHF] = max(peaks(tempidxHF));
    f_HF = F(idxpeaks(tempidxHF(idxHF)));
    
    if isempty(f_LF)
        intLF = find(F>=0.03 & F<=0.25);
        [peaks,idxpeaks] = findpeaks(PSD);
        tempidxLF = find(ismember(idxpeaks,intLF));
        [P_LF,idxLF] = max(peaks(tempidxLF));
        f_LF = F(idxpeaks(tempidxLF(idxLF)));
        
        intHF = find(F>0.25 & F<=0.50);
        tempidxHF = find(ismember(idxpeaks,intHF));
        [P_HF,idxHF] = max(peaks(tempidxHF));
        f_HF = F(idxpeaks(tempidxHF(idxHF)));
    end
    if ~isempty(f_LF) && ~isempty(f_HF)
        condition = 0;
    else
        optimalmodelorder = optimalmodelorder+1;
    end
end

% 2) variabiity indices
meanRR = mean(RR);
SDRR = std(RR);
rMSSDRR = sqrt(mean(diff(RR).^2));
pNN10 = numel(find(abs(diff(RR))>0.010))/length(diff(RR))*100;
pNN20 = numel(find(abs(diff(RR))>0.020))/length(diff(RR))*100;
pNN30 = numel(find(abs(diff(RR))>0.030))/length(diff(RR))*100;
pNN40 = numel(find(abs(diff(RR))>0.040))/length(diff(RR))*100;
pNN50 = numel(find(abs(diff(RR))>0.050))/length(diff(RR))*100;
pNN60 = numel(find(abs(diff(RR))>0.060))/length(diff(RR))*100;
pNN70 = numel(find(abs(diff(RR))>0.070))/length(diff(RR))*100;
pNN80 = numel(find(abs(diff(RR))>0.080))/length(diff(RR))*100;
pNN90 = numel(find(abs(diff(RR))>0.090))/length(diff(RR))*100;
pNN100 = numel(find(abs(diff(RR))>0.1))/length(diff(RR))*100;

% 3) irregularity indices
%{
saen1 = SampEn2(1,0.25*std(RR),RR);
if saen1 == Inf
    saen1 = SampEn2(1,0.35*std(RR),RR);
end
saen2 = SampEn2(2,0.25*std(RR),RR);
if saen2 == Inf
    saen2 = SampEn2(2,0.35*std(RR),RR);
end
%}

% 4) symbolic analysis
[V0, V1, LV2, UV2] = SymbAnal(RR,6,1);

%% P-QRS-T based measures

%%

% 1) complexity of QRS-T complex via PCA and shape analysis
% Tend = round(min(diff(QRS))*0.8);
RRmin = min(diff(QRS));
% QRSTMatrix = zeros(length(QRS)-2,Tend+0.1*fs+1);
QRSTMatrix = zeros(length(QRS)-2,round(0.3*RRmin)+round(0.7*RRmin)+1);
for i = 2:length(QRS)-1
%     QRSTMatrix(i-1,:) = ecg(QRS(i)-0.1*fs:QRS(i)+Tend);
    QRSTMatrix(i-1,:) = ecg(QRS(i)-round(0.3*RRmin):QRS(i)+round(0.7*RRmin));
end

% complexity
QRSTMatrix0 = bsxfun(@minus, QRSTMatrix, mean(QRSTMatrix,2));
[U, S, V] = svd(QRSTMatrix0,'econ');
N = size(QRSTMatrix0,2);
A = U*S/sqrt(N);
Z = sqrt(N)*V';
cumvar = cumsum(S.^2)/sum(S.^2);
numbercomponents95percQRST = find(cumvar>0.95,1,'first');

% plot(Z(1:5,:)')
% fa = (0:size(Z,2)-1)*fs/size(Z,2);
% [~,posindexmax] = max(abs(fft(Z')));
% componentwithdominatfrequency4to12Hz = find(fa(posindexmax) >=5 & fa(posindexmax) <= 12);
% QRSTMatrix0412 = A(:,componentwithdominatfrequency4to12Hz)*pinv(A(:,componentwithdominatfrequency4to12Hz))*QRSTMatrix0;
% QRSTMatrix412 = bsxfun(@plus, QRSTMatrix0412, mean(QRSTMatrix,2));
% plot(QRSTMatrix412')
% pause
% close
% for i = 1:5
% % [pa,fa]=pburg(Z(i,:),20,[],fs,'onesided');
% pa = abs(fft(Z(i,:)));
% plot(fa,pa)
% xlim([0 12])
% pause
% hold on
% end

template0 = A(:,numbercomponents95percQRST)*Z(numbercomponents95percQRST,:);
QRSTtemplate = bsxfun(@plus, template0, mean(QRSTMatrix,2));

% shape analysis
count = 0;
Dist = ones(size(QRSTMatrix,1)-1,size(QRSTMatrix,1))*(-2);
for i = 1:size(QRSTMatrix,1)-1
    for j = i+1:size(QRSTMatrix,1)
        Dist(i,j) = acos(QRSTMatrix(i,:)/norm(QRSTMatrix(i,:))*QRSTMatrix(j,:)'/norm(QRSTMatrix(j,:)));
        count = count+1;
    end
end

SimQRST1 = numel(find(Dist>-2 & Dist<pi/2))/count;
SimQRST2 = numel(find(Dist>-2 & Dist<pi/4))/count;
SimQRST3 = numel(find(Dist>-2 & Dist<pi/8))/count;
SimQRST4 = numel(find(Dist>-2 & Dist<pi/16))/count;
SimQRST5 = numel(find(Dist>-2 & Dist<pi/32))/count;

% 2) complexity of P-wave via PCA and shape analysis
PMatrix = zeros(length(QRS)-2,round(0.17*fs+1));
for i = 2:length(QRS)-1
    PMatrix(i-1,:) = ecg(QRS(i)-0.22*fs:QRS(i)-0.05*fs);
end

% complexity 
PMatrix0 = bsxfun(@minus, PMatrix, mean(PMatrix,2));
% [U, S, V] = svd(PMatrix0,'econ');
[~, S, ~] = svd(PMatrix0,'econ');
% N = size(PMatrix0,2);
% A = U*S/sqrt(N);
% Z = sqrt(N)*V';
cumvar = cumsum(S.^2)/sum(S.^2);
numbercomponents95percP = find(cumvar>0.95,1,'first');

% shape analysis
count = 0;
Dist = ones(size(PMatrix,1)-1,size(PMatrix,1))*(-2);
for i = 1:size(PMatrix,1)-1
    for j = i+1:size(PMatrix,1)
        Dist(i,j) = acos(PMatrix(i,:)/norm(PMatrix(i,:))*PMatrix(j,:)'/norm(PMatrix(j,:)));
        count = count+1;
    end
end

SimP1 = numel(find(Dist>-2 & Dist<pi/2))/count;
SimP2 = numel(find(Dist>-2 & Dist<pi/4))/count;
SimP3 = numel(find(Dist>-2 & Dist<pi/8))/count;
SimP4 = numel(find(Dist>-2 & Dist<pi/16))/count;
SimP5 = numel(find(Dist>-2 & Dist<pi/32))/count;


% To-DO:
% work also on dRR series
% symbolic analysis?
% P-wave absence


%% AA-based measures

%% single-lead AA extraction
% AA = ecg;
% for i = 2:length(QRS)-1
%     AA(QRS(i)-round(0.3*RRmin):QRS(i)+round(0.7*RRmin)) = ecg(QRS(i)-round(0.3*RRmin):QRS(i)+round(0.7*RRmin))-QRSTtemplate(i-1,:)';
% end

% AA(QRS(end)-round(0.3*RRmin):end) = [];
% AA(1:QRS(1)+round(0.7*RRmin)) = [];
% 
% plot([ecg(QRS(1)+round(0.7*RRmin)+1:QRS(end)-round(0.3*RRmin)-1) AA])
% pause
% close

% build atrial activity signal
AA = [];
for i = 2:length(QRS)-1
    AA = [AA;QRSTMatrix(i-1,:)-QRSTtemplate(i-1,:)];
end

% AA = [];
% for i = 2:length(QRS)-1
%     AA = [AA;QRSTMatrix412(i-1,:)];
% end

AA = AA';

% correct for misalignments
correctedAA=[];
for k = 2:length(QRS)-2
    c1=AA(end,k-1);
    c2=AA(1,k);
    N=length(QRS(k-1)+round(0.7*RRmin):QRS(k)-round(0.3*RRmin));
    bet=(ecg(QRS(k)-round(0.3*RRmin))-ecg(QRS(k-1)+round(0.7*RRmin))-c2+c1)/(N-1);
    alp=ecg(QRS(k-1)+round(0.7*RRmin))-c1-bet;
    vec=ecg(QRS(k-1)+round(0.7*RRmin):QRS(k)-round(0.3*RRmin))-(alp+bet*[1:N])';
    correctedAA=[correctedAA;vec(2:end-1);AA(:,k)];
end

% correctedAA=[];
% for k = 2:length(QRS)-2
%     correctedAA=[correctedAA;ecg(QRS(k-1)+round(0.7*RRmin)+1:QRS(k)-round(0.3*RRmin)-1);AA(:,k)];
% end

% smooth QRS residues
correctedAA = [zeros(length(1:QRS(1)+round(0.7*RRmin))-1,1);correctedAA;zeros(length(QRS(end-2)+round(0.7*RRmin)+1:length(ecg)),1)];
if min(diff(QRS))> round(0.1*fs)+9
    [Xar] = LSAR_original(correctedAA',QRS,fs);
    clear correctedAA
    correctedAA = Xar';
end
correctedAA(QRS(end-2)+round(0.7*RRmin)+1:length(ecg)-1) = [];
correctedAA(1:QRS(1)+round(0.7*RRmin)-1) = [];


% highpass filter to reduce T-wave residues
[b, a] = cheby2(3,20,3/(fs/2),'high');
filteredAA = filtfilt(b,a,correctedAA);

% compute spectral concentration and dominant frequency
[SCN, DF] = SC(filteredAA,fs);

% f-wave power, rhe, sample entropy, based on both filteredAA and main
% atrial wave signals
[fWP, fWPMAW, rhe,  SampEn, SampEnMAW, MAW] = SampEnIndicesAAsignals_NewVersion(filteredAA,fs);

% organization index based on 3 main peaks in 4-12 Hz
OI = OrganizationIndex(filteredAA,fs,3);

% spectral entropy
SE = SpectralEntropy(filteredAA);

% f-wave amplitude (findpeaks+cluster?)
[peakFWA] = findpeaks(MAW);

% outlier detection
Q1 = prctile(peakFWA,25);
Q3 = prctile(peakFWA,75);
IQR = Q3-Q1;
upperbound = Q3+1.5*IQR;
lowerbound = Q1-1.5*IQR;
% l(peakFWA>upperbound) = [];
peakFWA(peakFWA>upperbound) = [];
% l(peakFWA<lowerbound) = [];
peakFWA(peakFWA<lowerbound) = [];
% x = 1:length(MAW);
% plot(x,MAW,x(l),peakFWA,'or')
FWA = median(peakFWA);


%% create parameter vector
%{
parameters=[P_LF, f_LF, P_HF, f_HF, meanRR, SDRR, rMSSDRR, pNN10, pNN20, ...
        pNN30, pNN40, pNN50, pNN60, pNN70, pNN80,pNN90, pNN100,  ...
        V0, V1, LV2, UV2, numbercomponents95percQRST, SimQRST1, SimQRST2, ...
        SimQRST3, SimQRST4, SimQRST5, numbercomponents95percP, SimP1, SimP2, ...
        SimP3, SimP4, SimP5, SCN, DF, fWP, fWPMAW, rhe,  SampEn, SampEnMAW, ...
        OI, SE, FWA, differenceMediaRRmeanRR];
%}
clear parameters

parameters(1)=P_LF(1);
parameters(2)=f_LF(1); 
parameters(3)=P_HF(1);
parameters(4)=f_HF(1);
parameters(5)=meanRR(1);
parameters(6)=SDRR(1);
parameters(7)=rMSSDRR(1);
parameters(8)=pNN10(1);
parameters(9)=pNN20(1);
parameters(10)=pNN30(1);
parameters(11)=pNN40(1);
parameters(12)=pNN50(1);
parameters(13)=pNN60(1);
parameters(14)=pNN70(1);
parameters(15)=pNN80(1);
parameters(16)=pNN90(1);
parameters(17)=pNN100(1);
parameters(18)=V0(1);
parameters(19)=V1(1);
parameters(20)=LV2(1);
parameters(21)=UV2(1);
parameters(22)=numbercomponents95percQRST(1);
parameters(23)=SimQRST1(1);
parameters(24)=SimQRST2(1);
parameters(25)=SimQRST3(1);
parameters(26)=SimQRST4(1);
parameters(27)=SimQRST5(1);
parameters(28)=numbercomponents95percP(1);
parameters(29)=SimP1(1);
parameters(30)=SimP2(1);
parameters(31)=SimP3(1);
parameters(32)=SimP4(1);
parameters(33)=SimP5(1);
parameters(34)=SCN(1);
parameters(35)=DF(1);
parameters(36)=fWP(1);
parameters(37)=fWPMAW(1);
parameters(38)=rhe(1);
parameters(39)=SampEn(1);
parameters(40)=SampEnMAW(1);
parameters(41)=OI(1);
parameters(42)=SE(1);
parameters(43)=FWA(1);
parameters(44)=differenceMediaRRmeanRR(1);

