function TWARes = TWAbyMMAOnAFile
% TWAbyMMAOnAFile.m
% Author: Alexander Khaustov; alexander dot khaustov at gmail dot com 
% Copyright (C) 2008 St.-Petersburg Institute of Cardiological Technics (Incart), www.incart.ru
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).

    clear global TWARes Align CurrAvg
    global Param freq ecg q s stlen TWARes Align CurrAvg
    
    Param.Interval = 15;    % seconds
    intvLength = Param.Interval * freq; 

    TWARes.successfull = false;        
    TWARes.max_lead = 1;  
    TWARes.max_index = 1;
    
    [CurrAvg, Align] = InitializeMMAavg(ecg, q, s, stlen);
    if (isempty(CurrAvg))
        return;
    end;
    
    nextinterval = intvLength;
    StartQRSInd = 1;
    leads = size(ecg, 2);
    NInvalid = zeros(1, leads);

    for i = 1:length(q)
        ev_odd = mod(i, 2) + 1;
            
        if i > length(Align.fid)
            Align = AlignSingleBeat(ecg, q(i), s(i), Align, CurrAvg(ev_odd, Align.lead));
        end;
        
        for lead = 1:leads
            if Align.valid(i, lead)
                Align.amp(i, lead) = mean(ecg(Align.fidQRS(i) - Align.q2f : Align.fidQRS(i) + Align.f2s, lead));  % align to zero mean on QS segment
    
                [CurrAvg(ev_odd, lead).qs_avg, invalidQS] = MMAadd(CurrAvg(ev_odd, lead).qs_avg, ...
                    ecg(Align.fidQRS(i) - Align.q2f : Align.fidQRS(i) + Align.f2s, lead) - Align.amp(i, lead));
                
                [CurrAvg(ev_odd, lead).st_avg, invalidST] = MMAadd(CurrAvg(ev_odd, lead).st_avg, ...
                    ecg(Align.fid(i) + Align.f2s : Align.fid(i) + Align.f2s + Align.st, lead) - Align.amp(i, lead));                
                
                if (invalidQS > 0.1 * (Align.q2f + Align.f2s + 1) || invalidST > 0.1 * Align.st)
                    NInvalid(lead) = NInvalid(lead) + 1;
                end;
            else
                NInvalid(lead) = NInvalid(lead) + 1;
            end;
        end;
        
        if (i >= length(q) || q(i + 1) > nextinterval)
            intv = nextinterval / intvLength;
            
            for lead = 1:leads
                TWARes.res(intv).validleads(lead) = (NInvalid(lead) < 0.1 * (i - 1 - StartQRSInd));
            end;                                    
            TWARes.res(intv).StartQRSInd =  StartQRSInd;
            TWARes.res(intv).EndQRSInd = i;
            
            AddMMARes(intv);
          
            nextinterval = nextinterval + intvLength;
            StartQRSInd = i + 1;
            NInvalid = zeros(1, leads);
        end;
        
    end;
    
    if (TWARes.successfull)
        [m mi] = max(TWARes.VAltTrend);
        [TWARes.maxVAlt TWARes.max_lead] = max(m);
        TWARes.max_index = mi(TWARes.max_lead);
    end;

return;
    
function AddMMARes(intv)

global TWARes Align CurrAvg
leads = size(CurrAvg, 2);

TWARes.res(intv).Avg = CurrAvg;

for lead = 1:leads
    if (TWARes.res(intv).validleads(lead))
        df = TWARes.res(intv).Avg(1, lead).st_avg - TWARes.res(intv).Avg(2, lead).st_avg;
        [ma TWARes.res(intv).VAltPt(lead)] = max(abs(df));
%         s = std(df);
%           if (ma > 3 * s)
            TWARes.res(intv).VAlt(lead) = abs(mean(df)); %ma;
            TWARes.successfull = true;
%         else
%             TWARes.res(intv).VAlt(lead) = 0;
%         end;
    else
        TWARes.res(intv).VAlt(lead) = 0;
    end;
end;

TWARes.VAltTrend(intv, :) = TWARes.res(intv).VAlt;

return;