function lpcpoly = lsf2poly(lsf)
% Convert a vector of line spectral frequencies to the equivalent linear
% prediction polynomial.
%
% If A(z) = 1 + a1*z^-1 + ... + am*z^-M is
% a linear prediction polynomial, the corresponding line
% spectral frequencies are the angles of the roots in the upper half
% plane of the two polynomials
% P_sym(z) = A(z) + (z^-(M+1))*A(z^-1) and
% P_asym(z) = A(z) - (z^-(M+1))*A(z^-1), ignoring the roots 1 and -1,
% which are always present and so don't need to be listed.
% P_sym and P_asym are known as the symmetric and asymmetric polynomials,
% respectively, and it is easy to see that P_sym(z) + P_asym(z) = 2*A(z).
% It can be shown that the roots of P_sym and P_asym all lie on the unit
% circle, so they are completely characterized by their angles.
% Furthermore, the roots of P_sym and P_asym interleave, which is to say
% that there is a root of P_sym between any pair of roots of P_asym, and
% vice versa.
% It is clear by inspection that 1 (with angle 0 in the complex plane) is
% a root of P_asym. Thus if we take the lsfs together with the angles
% of the conjugate roots in the lower half plane, plus 0 and pi, and put
% them in ascending order, the first (i.e., 0) and every second one after
% that will be the angles of roots of P_asym, and the remainder will be
% angles of roots of P_sym. We can then reconstruct P_sym and P_asym from
% their roots and obtain A by adding them and dividing by 2.
% P_asym always has a leading coefficient of -1, but the octave
% poly function always returns a result with a positive leading
% coefficent, so we negate the polynomial returned by the poly function.
% We discard the leading coefficient of the sum because the leading 1 and
% -1 coefficients of P_sym and P_asym always sum to zero.
% Finally we reverse the order of the lpc coefficients because they are
% traditionally presented with the constant term first, while the output
% of the poly function puts the constant term last.
% For compatibility with matlab, if the input is a matrix whose columns
% are lsfs, the output is a matrix whose rows are the corresponding lpc
% polynomials. If the input is either a row or column vector, the output
% is a row vector.
isreal(lsf) || error("lsfs have to be real");
(any(lsf < 0) || any(lsf > pi)) && error("lsfs lie between 0 and pi");
if(size(lsf,1) == 1)
lsf = lsf';
end
ncols = size(lsf,2);
lpcpoly = zeros(ncols,size(lsf,1)+1);
for j=1:ncols
rts = exp(i*lsf(:,j));
rts = [1; rts; -1; flipud(conj(rts))];
p_asym = real(poly(rts(1:2:end)));
p_sym = real(poly(rts(2:2:end)));
% For some reason octave makes the coefficients in poly([1 r conj(r)])
% be complex, with 0.0i imaginary part, although the coefficients of
% poly([r conj(r)]) are given as real.
lpcpoly(j,:) = fliplr((p_sym - p_asym)(2:end))/2;
end
end