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
