function [TWARes, Align] = TWASpectral(ecg, Align)
% TWASpectral.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).
% 
% TWA by Spectral Method on a given interval
% INPUT:
% ecg   ecg signal (could be multiple channels)
% Align alignment structure produced by AlignBeats
% Note the following definitions:
%       num_of_beats is the length of the analysis window, typically 128 beats
%       num_of_leads is the number of ECG leads
%       num_of_timepoints is the length of the Q-Tend complex
%
% OUTPUT:
% TWARes: structure containing the fields:
%       alt_series     TWA sequence array
%       avg_even (num_of_timepoints x num_of_leads): time average of even beats
%       avg_odd  (num_of_timepoints x num_of_leads): time average of odd beats
%       psd      (num_of_beats/2+1 x num_of_leads x num_of_timepoints): psd of timepoints within Q-Tend
%       lomb     (num_of_beats/2+1 x num_of_leads x num_of_timepoints): lomb psd of timepoints within Q-Tend
%       avg_psd  (num_of_beats/2+1 x num_of_leads): average psd across all timepoints in Q-Tend
%       significant    (1 x num_of_leads): 1 if the alternans value is statistically significant against the Param.RatioThreshold in the lead and 0 otherwise
%       VAlt     (1 x num_of_leads):
%       Ratio    (1 x num_of_leads):
%
% ecg: num_of_samples x number_of_leads


global Param

TWARes = [];

if ~isempty(strfind(Param.MethodForEctopy, 'lomb'))
    TWARes.lomb.successfull = false;
end;
if ~isempty(strfind(Param.MethodForEctopy, 'replace'))
    TWARes.replace.successfull = false;
end;
if ~isempty(strfind(Param.MethodForEctopy, 'differences'))
    TWARes.differences.successfull = false;
end;
    
for lead = 1:size(ecg, 2)

    if (~Align.validleads(lead))
        continue;
    end;
    
    if ~isempty(strfind(Param.MethodForEctopy, 'lomb'))
        [TWARes.lomb.at_lead(lead).series(:, :), TWARes.lomb.at_lead(lead).times] = ...
            CalcAltSeriesForLomb(ecg(:, lead), Align.fid + Align.f2s, Align.amp(:, lead), Align.st, Align.valid(:, lead));
    
        global beats
        len = beats;
        f = [0:len / 2] / len;
        for timept = 1:size(TWARes.lomb.at_lead(lead).series, 2)
            a = detrend(TWARes.lomb.at_lead(lead).series(:, timept));
            [b, prob] = lomb(TWARes.lomb.at_lead(lead).times', a, f);
            TWARes.lomb.psd(:, lead, timept) = b * sum(a .* a) / sum(b) / length(a);    % normalization to equal energy with simply psd
        end;
    end;
    
    names = cellstr({'replace', 'differences'});
    for i = 1:length(names)
        if ~isempty(strfind(Param.MethodForEctopy, names{i}))
            [TWARes.(names{i}).at_lead(lead).series(:, :), ae, ao] = ... 
                CalcAltSeries(ecg(:, lead), Align.fid + Align.f2s, Align.amp(:, lead), Align.st, Align.valid(:, lead), i == 2);
            if ~isempty(ae)
                TWARes.(names{i}).avg_even(:, lead) = ae;
                TWARes.(names{i}).avg_odd(:, lead) = ao;
            end;
            
            [TWARes.(names{i}).psd(:, lead, :)] = CalcPSD(squeeze(TWARes.(names{i}).at_lead(lead).series(:, :)), i == 2); 
        end;
    end;

    names = cellstr({'lomb', 'replace', 'differences'});
    for i = 1:length(names)
        if ~isfield(TWARes, names{i})
            continue;
        end;
        
        TWARes.(names{i}).avg_psd(:, lead) = CalcAvgPSD(squeeze(TWARes.(names{i}).psd(:, lead, :)));
        [TWARes.(names{i}).significant(lead) TWARes.(names{i}).VAlt(lead) TWARes.(names{i}).Ratio(lead)] = CalcValues(TWARes.(names{i}).avg_psd(:, lead));
        
    
        if (TWARes.(names{i}).significant(lead) && TWARes.(names{i}).Ratio(lead) > 3)
            TWARes.(names{i}).successfull = true;
        end;
    end;
    
end;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%