How to make that process faster? Well, firstly if you feel like having the time you could check the implementation, which I'm doing just right now.
(Warning: I need to truncate links as stackexchange does not let me have more than 2 links)
Code is here:
assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum
Giving a look at the package it already contains some .c code to fasten up some operations. That may ease the implementation in C of some part to later on import it via Python if needed.
I did
svn checkout https://subversion.assembla.com/svn/PySpectrum/
And I got:
~/spectrumlib/PySpectrum/trunk/src/cpp$ ls
init.py mydpss.c
mydpss.c implements some subroutines. Something could be added there maybe.
But let's take a look at the code that is actually being executed:
assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum/yulewalker.py
from correlation import CORRELATION
from levinson import LEVINSON
r = CORRELATION(X, maxlags=order, norm=norm)
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
return A, P, k
This is what is being executed. Let's check what is the slow part here.
So I modified yulewalker.py like:
import time
print "Correlation calculation"
t0 = time.time()
r = CORRELATION(X, maxlags=order, norm=norm)
print "Time elapsed: ", time.time()-t0
print "Levinson calculation"
t1 = time.time()
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
print "Time elapsed: ", time.time()-t1
return A, P, k
I ran your piece of code of stackexchange and I got:
./aryulespeed.py
Correlation calculation
Time elapsed: 1.50509214401
Levinson calculation
Time elapsed: 0.000110864639282
Time elapsed: 1.50538802147
[-2.7616739 3.81292317 -2.65578205 0.9248034 ]
So the correlation calculation is the slow part! aha!
This guy:
assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum/correlation.py
seems to be the fat guy.
Oh, look, a funny comment in the code:
Provides two correlation functions. :func:`CORRELATION` is slower than
:func:`xcorr`. However, the output is as expected by some other functions.
Ultimately, it should be replaced by :func:`xcorr`.
For real data, the behaviour of the 2 functions is identical. However, for
complex data, xcorr returns a 2-sides correlation.
Maybe we should try to interchange this CORRELATION for xcorr!
So I tried:
#! /usr/bin/env python
from pylab import *
import scipy.signal
from spectrum import *
import time
# The original imports
from spectrum.correlation import CORRELATION
from spectrum.levinson import LEVINSON
# The faster corr
from spectrum.correlation import xcorr
def superawesomearyule(X, order, norm='biased', allow_singularity=True):
assert norm in ['biased', 'unbiased']
import time
print "Correlation calculation in superawesomearyule"
t0 = time.time()
#r = CORRELATION(X, maxlags=order, norm=norm) # Commenting you, slowwww code!
r = xcorr(X, maxlags=order, norm=norm)
print "Time elapsed: ", time.time()-t0
print "Levinson calculation"
t1 = time.time()
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
print "Time elapsed: ", time.time()-t1
return A, P, k
# Create an AR model
a = [1, -2.7607, 3.8106, -2.6535, 0.9238]
# create some data based on these AR parameters
X = randn(1, 256000)
y = scipy.signal.lfilter([1], a, X)
# now, let us try to estimate the original AR parameters
t0 = time.time()
#AR, P, k = aryule(y[0], 4)
AR, P, k = superawesomearyule(y[0], 4)
print "Time elapsed: ", time.time()-t0
#%compare the results in ar_coeffs to the vector A.
print(AR)
And executed:
./aryulespeed.py
Correlation calculation in superawesomearyule
Time elapsed: 69.9396679401
Levinson calculation
Traceback (most recent call last):
File "./aryulespeed.py", line 38, in <module>
AR, P, k = superawesomearyule(y[0], 4)
File "./aryulespeed.py", line 24, in superawesomearyule
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
File "/usr/local/lib/python2.7/dist-packages/spectrum/levinson.py", line 112, in LEVINSON
if P <= 0 and allow_singularity==False:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
And...... oh noes :( it took 69 seconds!! and also crashed after that on levinson.
I think giving it a bit of love something may be done here. Maybe sending a mail to it's developer Thomas Cokelaer would help.
And I'm sorry but I can't invest more time in this issue right now. I hope this helps you in some way :)
I also found this github repo:
https://github.com/RhysU/ar
Which has an ar implementation with some python bindings. It may also help you.
Good luck!