function [TWARes] = TWA_MMA(Param, freq, ecg, q,r,s, stlen)
%OVERVIEW, This function checks for arrhythmia in the current analysis window across all
% leads. If no arrhythmia is detected it then calls sub functions which first align the beats in the analysis 
% window and find correlations of each beat with a template, to determine if the ecg data is clean enough for analysis.
% If alignment is successful, TWAs are evaluated using the MMA method
% coupled with a surrogate non parametric reshuffling method used for
% reducing false positive TWA detections.
%
% INPUTS        MANDATORY           DESCRIPTION
%               Param               Structure contains parameters used for
%                                   computing TWAs.
%
%               freq                Sampling frequency of the ecg data,
%                                   needs to be 1000 Hz.
%
%               ecg                 N by M array of ecg data. Contains M
%                                   channels of ECG with N datapoints each.
%
%               q                   array of q point locations
%
%               r                   array of r point locations
%
%               s                   array of s point locations
%
%               stlen               scalar, estimate of s-t segment length
%
% OUTPUT
%
%               TWARes              Structure contains the result for TWA
%                                   computation. Contains the following fields,
%
%               .VAlt               TWA amplitude estimates for the analysis
%                                   window.
%
%               .VAlt_Sig           TWA amplitude estimates which are statistically
%                                   significant compared to the noise threshold
%                                   estimate.
%
%               .noise_median       median of the gamma distribution used to model the noise.
%
%               .noise_95           95th percentile of the gamma distribution
%                                   used to model the noise.
%
%               .VAltPt             The location of the point with maximum
%                                   difference between the average even and odd beats.
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Shamim Nemati, 
%       editted by Ismail Sadiq on 10/26/2019. 
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.  

TWARes.successfull = false;

WIN_SIZE = Param.Interval; %e.g., 60beats
OVERLAP_SIZE = floor(WIN_SIZE/2);
TWARes=[];
Num_Ch = size(ecg, 2);

Align.fidBase = []; Align.fid = []; Align.valid = [];

if (length(r) < WIN_SIZE) % need at least one window of beats 
    disp('TWAbyMMAOnAFile: Not enough beats...')
    TWARes.HR=NaN;
    TWARes.VAlt=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig=NaN*ones(1,Num_Ch);
    TWARes.VAltPt=NaN*ones(1,Num_Ch);
    TWARes.Noise_Median= NaN*ones(1,Num_Ch);TWARes.Noise_95 = NaN*ones(1,Num_Ch);
    return;
end

i = 0;
for j = 1:(WIN_SIZE-OVERLAP_SIZE):length(r)-WIN_SIZE+1
    i=i+1;
    disp(['Analysis window index = ' num2str(i)])
    beats = j:j+WIN_SIZE-1;
    
    % Check for arrhythmia
    
    cur_r = r(beats); % Check current window of beats for arrhythmia 
    features = AF_features(diff(cur_r),freq);
    Index_test = SVM_AFdetection_withoutTrainingModel(features,1);
    if (Index_test)
        disp('Arrhythmia detected...')
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :) = NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        % If arrhythmia detected continue to next window
        continue;
    end
    
    
    try
        % Align beats and check validity for analysis before computing TWAs
        [Align hrt_rate] = Align_Beats(beats, q(beats), r(beats), s(beats), stlen, ecg, 1000, Param);
    catch
        % If failed to align due to catch continue to next window.
        disp('TWAbyMMAOnAFile: Catch says alignment failed...')
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch);TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :) = NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    end
    if (isempty(Align.fid))
        disp('TWAbyMMAOnAFile: alignment failed ...');
        TWARes.HR(i)=NaN;
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch); TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :)=NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    else
        disp('TWAbyMMAOnAFile: alignment succeeded ...');
    end
    
    disp('TWAbyMMAOnAFile: looking for alternans ...');
    TWARes.HR(i) = hrt_rate;
    try
        [TWARes.VAlt(i, :), TWARes.VAlt_Sig(i, :), TWARes.Noise_Median(i,:),TWARes.Noise_95(i,:),...
            TWARes.VAltPt(i,:)]= TWA_by_MMA(Param, TWARes, ecg, Align); % *********************** Estimates T-wave alternans from the beat Matrix
    catch
        disp('TWAbyMMAOnAFile: search for alternans failed ')
        TWARes.VAlt(i, :)=NaN*ones(1,Num_Ch); TWARes.VAlt_Sig(i, :)=NaN*ones(1,Num_Ch);
        TWARes.VAltPt(i,:)=NaN*ones(1,Num_Ch);
        TWARes.Noise_Median(i, :)=NaN*ones(1,Num_Ch);TWARes.Noise_95(i, :) = NaN*ones(1,Num_Ch);
        continue;
    end
end


TWARes.successfull = true;

return;
