Search code examples
pythonsignal-processing

spectrum examples don't work?


I downloaded the marple data from here:

http://www.ece.rice.edu/dsp/courses/elec532/DATA/

I then tried to run examples from here:

http://thomas-cokelaer.info/software/spectrum/html/user/ref_param.html#module-burg

At first, I thought maybe the data was in the wrong format, so I converted the data to remove any leading whitespace and convert the space in between to a single comma with:

cat marple.ascii | sed 's/^[ \t]*//;s/[ \t]*$//' > output.txt
sed 's/ \{1,\}/,/g' output.txt > marple_comma.dat

and then reading the file with:

marple_data = numpy.loadtxt('marple_comma.dat', delimiter = ',', dtype = numpy.complex128)

but I get a different error:

TypeError: only length-1 arrays can be converted to Python scalars

I tried the example below. But I get this error. What is wrong?

runfile('/home/idf/marple.py', wdir='/home/idf')
Traceback (most recent call last):

  File "<ipython-input-14-e02af4ddfead>", line 1, in <module>
    runfile('/home/idf/marple.py', wdir='/home/idf')

  File "/home/idf/anaconda/lib/python2.7/site-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 682, in runfile
    execfile(filename, namespace)

  File "/home/idf/anaconda/lib/python2.7/site-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 78, in execfile
    builtins.execfile(filename, *where)

  File "/home/idf/marple.py", line 17, in <module>
    p()

  File "/home/idf/anaconda/lib/python2.7/site-packages/spectrum/arma.py", line 314, in __call__
    ma_params, rho = ma(self.data, self.ma_order, self.ar_order)

  File "/home/idf/anaconda/lib/python2.7/site-packages/spectrum/arma.py", line 377, in ma
    a, rho, _c = yulewalker.aryule(X, M, 'biased')   #! Eq. (10.5)

  File "/home/idf/anaconda/lib/python2.7/site-packages/spectrum/yulewalker.py", line 110, in aryule
    r = CORRELATION(X, maxlags=order, norm=norm)

  File "/home/idf/anaconda/lib/python2.7/site-packages/spectrum/correlation.py", line 131, in CORRELATION
    r[k-1] = sum / float(N)

ValueError: setting an array element with a sequence.





import numpy as np

from pylab import *
from spectrum import *

#marple_data = np.genfromtxt('marple.ascii')
marple_data = numpy.loadtxt('marple.ascii') 

p = pma(marple_data, 15, 30, NFFT=4096)
p()
p.plot(sides='centerdc')

Solution

  • I think your input file contains the real and imaginary parts as space-delimited numbers on each row: marple.ascii:

      1.349839091    2.011167288
     -2.117270231    0.817693591
    ...
     -0.895898521   -0.364855707
    

    If you use loadtxt to read in the whole array it isn't clever enough to realise that the two columns are the real and imaginary parts and form a single complex number from each row.

    If this is the case, you can read it in as a 64x2 array of floats (ascii, below) and then build your complex array, data from it:

    In [1]: import numpy as np
    In [2]: ascii = np.loadtxt('marple.ascii')
    In [3]: data = ascii[:,0] + 1j*ascii[:,1]
    In [4]: data
    array([ 1.34983909+2.01116729j, -2.11727023+0.81769359j,
           -1.78642166-1.29169893j,  1.16223633-1.48259807j,
           ...
           -0.76273143+0.40897125j, -0.89589852-0.36485571j])
    

    An alternative to line [3] if you have a lot of data and don't want to create a whole new array is:

    data = ascii.view(dtype=np.complex128)[:,0]