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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%