% 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) 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);