频谱示例不起作用?

spectrum examples don't work?

我从这里下载了 marple 数据:

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

然后我尝试 运行 来自这里的示例:

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

起初,我认为可能是数据格式错误,所以我转换数据以删除任何前导白色space并将中间的space转换为单个逗号:

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

我认为您的输入文件包含每行以 space 分隔的数字形式的实部和虚部: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]