频谱示例不起作用？

Question

我从这里下载了 marple 数据：

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

然后我尝试从这里运行示例：

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

起初，我认为可能是数据格式错误，因此我将数据转换为删除任何前导空格并将其间的空格转换为单个逗号：

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

然后读取文件：

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

但我得到一个不同的错误：

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

我尝试了下面的示例。但我得到这个错误。 怎么了？

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')

Answer 1

我认为您的输入文件包含实部和虚部，每行以空格分隔的数字：marple.ascii:

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

如果您使用loadtxt 读取整个数组，那么意识到两列是实部和虚部并从每一行形成一个复数是不够聪明的。

如果是这种情况，您可以将其读取为 64x2 浮点数组（ascii，如下），然后从中构建您的复杂数组 data：

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])

如果您有大量数据并且不想创建一个全新的数组，则第 [3] 行的替代方法是：

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