PBURG example ================ Here is another method to estimate an AR model, based on :func:~spectrum.burg.arburg . This example is inspired by an example found in Marple book. This is very similar to the previous example, where you will find more explanation (see yule-Walker tutorial). .. doctest:: from pylab import * import scipy.signal from spectrum import * # Define AR filter coefficients a = [1, -2.2137, 2.9403, -2.1697, 0.9606]; .. doctest:: [w,H] = scipy.signal.freqz(1, a, 256) Hp = plot(w/pi, 20*log10(2*abs(H)/(2.*pi)),'r') .. doctest:: x = scipy.signal.lfilter([1], a, randn(256)) AR, rho, ref = arburg(x, 4) .. doctest:: PSD = arma2psd(AR, rho=rho, NFFT=512) PSD = PSD[len(PSD):len(PSD)/2:-1] plot(linspace(0, 1, len(PSD)), 10*log10(abs(PSD)*2./(2.*pi))) xlabel('Normalized frequency (\times \pi rad/sample)') .. plot:: :width: 80% from pylab import * import scipy.signal from spectrum import * # Define AR filter coefficients a = [1, -2.2137, 2.9403, -2.1697, 0.9606]; [w,H] = scipy.signal.freqz(1, a, 256) Hp = plot(w/pi, 20*log10(2*abs(H)/(2.*pi)),'r') x = scipy.signal.lfilter([1], a, randn(256)) AR, rho, ref = arburg(x, 4) PSD = arma2psd(AR, rho=rho, NFFT=512) PSD = PSD[len(PSD):len(PSD)/2:-1] plot(linspace(0, 1, len(PSD)), 10*log10(abs(PSD)*2./(2.*pi))) xlabel('Normalized frequency (\times \pi rad/sample)') ylabel('One-sided PSD (dB/rad/sample)') legend(['PSD of model output','PSD estimate of x']) .. plot:: :width: 80% from pylab import * import scipy.signal from spectrum import * # Define AR filter coefficients a = [1, -2.2137, 2.9403, -2.1697, 0.9606]; [w,H] = scipy.signal.freqz(1, a, 256) Hp = plot(w/pi, 20*log10(2*abs(H)/(2.*pi)),'r') x = scipy.signal.lfilter([1], a, randn(256)) p = pburg(x, 4, sampling=2) p() p.plot()