function [ap,aq] = lsptd(wa,p)
% LSPTD   Line Spectrum Pair (LSP) decomposition in the Time Domain
%    [AP,AQ] = LSPTD(WA,P) Finds the coefficients for two optimal
%    predictors equivalent to the LSP-decompositions. The coefficients
%    are AP=[1 AP(2) ... AP(P+1)] and AQ=[1 AQ(2) ... AQ(P+1)].
%    Note: Trivial zeros are included.
%
%    This function is an appendix to the article "A time domain 
%    reformulation of linear prediction equivalent to the LSP 
%    decomposition", by Bäckström, Alku and Kleijn, presented
%    at ICASSP'02.
%
%    See also: LSP

% Tom Bäckström / 11.10.2001

% calculate autocorrelation

R = xcorr(wa(:),p+1);


% generate special autocorrelation matrix

Rm1 = toeplitz(2*R(1+p+(1:p)) + R(p+(1:p)) + R(2+p+(1:p)));
Rm2 = toeplitz(2*R(1+p+(1:p)) - R(p+(1:p)) - R(2+p+(1:p)));


% generate special autocorrelation vector

Rv1 = -2*(R(1+p+(1:p)) + R(2+p+(1:p)));
Rv2 = -2*(R(1+p+(1:p)) - R(2+p+(1:p)));


% invert and multiply

hv1 = Rm1 \ Rv1;
hv2 = Rm2 \ Rv2;


% generate transfer function

ap = [1 zeros(1,p)] + [hv1(1:p)' 0]/2 + [0 hv1(1:p)']/2;
aq = [1 zeros(1,p)] + [hv2(1:p)' 0]/2 - [0 hv2(1:p)']/2;


% scale a0 to unity

ap = ap / ap(1);
aq = aq / aq(1);


