from spectrum import arma_estimate, arma2psd, marple_data
from pylab import plot, axis, xlabel, ylabel, grid, log10
ar, ma, rho = arma_estimate(marple_data, 15, 15, 30)
psd = arma2psd(ar, ma, rho=rho, NFFT=4096)
plot(10*log10(psd/max(psd)))
axis([0, 4096, -80, 0])
xlabel('Frequency')
ylabel('power (dB)')
grid(True)