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