function [alt_s, avg_even, avg_odd] = CalcAltSeries(ecg, fid, amp, STLen, valid, calc_from_diff)
% CalcAltSeries.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).
%
% calculates the alternans series for different points within st segment.
% If a given beat is marked as invalid it's replaced by the average
% even/odd beat
% this is the routine used for Param.MethodForEctopy = 'replace' and
% Param.MethodForEctopy = 'differences'
%
% normally the series is formed from sequential amplitudes of equally timed
% points within ST segment, but if the parameter calc_from_diff is
% set - series is calculated as differences between amplitudes in adjacent
% beats which performs natural detrending of the data that is later subject
% to power spectral estimation and supresses low frequencies
%
% INPUT:
% see AlignBeats.m
% OUTPUT:
% alt_s (num_of_beats x num_of_leads x num_of_timepoints): TWA series
% avg_even (num_of_timepoints x num_of_leads)
% avg_odd (num_of_timepoints x num_of_leads)
avg_even = [];
avg_odd = [];
invalid_exist = ~isempty(find(valid==0));
if (invalid_exist)
vinds = find(valid);
odd = vinds(find(mod(vinds(:), 2))); % odd indices of valid QRS
even = vinds(find((1 - mod(vinds(:), 2)))); % even indices of valid QRS
end;
alt_s = zeros(length(fid) - 1, STLen);
for timept = 1:STLen
% whenever there are invalid complexes compute odd and even
% averages to replace those
if (invalid_exist)
avg_odd(timept) = mean(ecg(fid(odd) + timept) - amp(odd));
avg_even(timept) = mean(ecg(fid(even) + timept) - amp(even));
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:(length(fid) - 1)
if (valid(i))
first = ecg(fid(i) + timept) - amp(i);
elseif mod(i, 2)
first = avg_odd(timept);
else
first = avg_even(timept);
end;
if (~calc_from_diff)
alt_s(i, timept) = first;
else
if (valid(i + 1))
second = ecg(fid(i + 1) + timept) - amp(i + 1);
elseif mod(i + 1, 2)
second = avg_odd(timept);
else
second = avg_even(timept);
end;
alt_s(i, timept) = first - second;
end;
end;
end;
return;