% function r = specfac(q)
%
% Find the minimum phase spectral factor of the polynomial whose
% coefficients are q(n).
%
% Kevin Amaratunga
% 2 March, 1998
function r = specfac(q)
L = length(q);
% Check that q is admissible.
if rem(L,2) == 0
warning('q[n] must have odd length.')
end
for k = 1:(L+1)/2
if q(k) ~= q(L-k+1)
warning('q[n] must be symmetric.')
end
end
d = (L - 1) / 2; % Delay that needs to be applied
% in order to make q causal.
% Compute the DFT of q.
M = 512; % Size of the DFT (make this even.)
qp = zeros(1,M);
qp(1:L-d) = q(d+1:L);
qp(M-d+1:M) = q(1:d);
Q = fft(qp);
% Q must be real, assuming that q is symmetric. Check that Q is also positive.
if nnz(Q) < M
warning('Q(z) has zeros on the unit circle. Cannot take the logarithm of 0.')
end
if nnz(real(Q) < 0) > 1
warning('Negative Q(w) is a bad sign. q[n] cannot be an autocorrelation.')
end
% Compute the logarithm.
Qhat = log(Q);
% Compute the cepstrum.
qhat = ifft(Qhat);
% Find the causal part.
rhat = zeros(1,M);
rhat(1) = qhat(1) / 2;
rhat(2:M/2) = qhat(2:M/2);
% Determine the DFT of the causal part.
Rhat = fft(rhat);
% Take the exponent.
R = exp(Rhat);
% Take the inverse DFT to get the spectral factor.
r = ifft(R);
% Pick out the right number of coefficients.
N = (L + 1) / 2;
if abs(r(N+1:M)) > 1e-10
warning('Unexpected non-zero coefficients in r[n].')
end
r = r(1:N);
% Finally, r may have a small imaginary part due to roundoff error,
% so just retain the real part.
if (imag(r)) > 1e-10
warning('Unexpected imaginary part in r[n].')
end
r = real(r);