TWAnalyser - A T-wave Alternans Detector 1.0.0

function drawres = TWA(varargin)
% TWA.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).
% 
% The analysis itself. Reads the file, checks for annotations, performs
% basic checks, filtration and analysis

    global record
    if ~nargin
        [fname, pathname] = uigetfile('*.*');
        if ~fname
            drawres = false;
            return;
        end;
        record = fname(1:strfind(fname, '.') - 1);
        drawres = TWA(record);
        return;
    else
        record = varargin{1};
    end;    
    
    Annotator = 'vf';

    if (~exist(strcat(varargin{1}, '.', Annotator),'file'))      
        %   ecgpuwave analyses an ECG signal from the specified record, detecting the QRS complexes and locating the 
        %   beginning, peak, and end of the P, QRS, and ST-T waveforms. 
        Annotator = 'pu';
        if (~exist(strcat(varargin{1}, '.', Annotator),'file'))      
            ss = sprintf('ecgpuwave -r %s -a pu', varargin{1});
            disp(ss);
            system(ss);
        end
    end;
        
% modified by Shamim to include NSAMP and TSTART
    if (nargin > 2)
        NSAMP = varargin{2};
        TSTART = varargin{3};
        readsignalfromwfdb(varargin{1}, Annotator, NSAMP,TSTART);
    else
        readsignalfromwfdb(varargin{1}, Annotator);
    end
    
    global ecg q s qs stend freq beats stlen Param

    %   removing last QRS to avoid problems with ecg length
    s = s(1:length(s) - 1);
    q = q(1:length(s));

    beats = 129;
    
% interval of fiducial point adjustment (radius): depends on possible Q, S
% detection errors and QT variation; seems to be even
% worse for MMA:
    if strcmp(Param.Metric, 'SM')
        Param.stAdjIntv = floor(30 * freq / 1000); %    see twa69, vcgamp1clean for the worst case; 
    else 
        Param.stAdjIntv = floor(40 * freq / 1000);
    end;

    Param.corrQRS = 0.96;   %   twa95
    Param.corrT = 0.8;
    
    Param.RatioThreshold = 3;   % for spectral method

    
    if (strcmp(Param.Metric, 'SM') && length(s) < 0.9 * beats) % twa43, twa84.. should be 0.9 * beats)
        disp('Not enough QRS');
        drawres = false;
        return;
    end;

    drawres = true;
    
    s = TWAUpdateArray(s, beats);
    q = TWAUpdateArray(q, beats);
    qs = TWAUpdateArray(qs, beats);

    stend = TWAUpdateArray(stend, beats);
    stlen = ApproximateSTLen(s, stend);

    %   pre-filtration of the signal
    Param.Filter = '';  % insert filter titles here
    ecg = FilterForTWA(ecg, freq);
    
    if (strcmp(Param.Metric, 'SM'))
        TWASpectralOnAFile;
    else
        TWAbyMMAOnAFile;
    end;

return;