Spectral analysis of a two frequencies signal¶
Context¶
Example¶
In the following example, we use most of the methods available to analyse an input signal made of the addition of two sinus and an additive gaussian noise
import numpy
import spectrum
from spectrum import tools
from numpy.testing import assert_array_almost_equal
import pylab
data = spectrum.marple_data
from pylab import *
nn = numpy.arange(200)
xx = cos(0.257*pi*nn) + sin(0.2*pi*nn) + 0.01*randn(size(nn));
def create_all_psd():
f = pylab.linspace(0, 1, 4096)
pylab.figure(figsize=(12,8))
# MA model
p = spectrum.pma(xx, 64,128); p(); p.plot()
"""
#ARMA 15 order
a, b, rho = spectrum.arma_estimate(data, 15,15, 30)
psd = spectrum.arma2psd(A=a,B=b, rho=rho)
newpsd = tools.cshift(psd, len(psd)//2) # switch positive and negative freq
pylab.plot(f, 10 * pylab.log10(newpsd/max(newpsd)), label='ARMA 15,15')
"""
# YULE WALKER
p = spectrum.pyule(xx, 7 , NFFT=4096, scale_by_freq=False); p.plot()
# equivalent to
# plot([x for x in p.frequencies()] , 10*log10(p.psd)); grid(True)
#burg method
p = spectrum.pburg(xx, 7, scale_by_freq=False); p.plot()
#pcovar
p = spectrum.pcovar(xx, 7, scale_by_freq=False); p.plot()
#pmodcovar
p = spectrum.pmodcovar(xx, 7, scale_by_freq=False); p.plot()
# correlogram
p = spectrum.pcorrelogram(xx, lag=60, NFFT=512, scale_by_freq=False); p.plot()
# minvar
p = spectrum.pminvar(xx, 7, NFFT=256, scale_by_freq=False); p.plot()
# pmusic
p = spectrum.pmusic(xx, 10,4, scale_by_freq=False); p.plot()
# pmusic
p = spectrum.pev(xx, 10, 4, scale_by_freq=False); p.plot()
# periodogram
p = spectrum.Periodogram(xx, scale_by_freq=False); p.plot()
#
legend( ["MA 32", "pyule 7", "pburg 7", "pcovar", "pmodcovar", "correlogram",
"minvar", "pmusic", "pev", "periodgram"])
pylab.ylim([-80,80])
create_all_psd()
Total running time of the script: (0 minutes 0.326 seconds)
Download Python source code:
plot_allpsd.py
Download IPython notebook:
plot_allpsd.ipynb