josef.pktd@gmail.com wrote:

On Sun, May 24, 2009 at 6:16 AM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:

...

Hi,

I have taken a look at the correlate function in scipy.signal. There are several problems with it. First, it is wrong on several accounts: - It assumes that the correlation of complex numbers corresponds to complex multiplication, but this is not the definition followed by most textbooks, at least as far as signal processing is concerned. - More significantly, it is wrong with respect to the ordering: it assumes that correlate(a, b) == correlate(b, a), which is not true in general.

I don't see this in the results. There was recently the report on the mailing list that np.correlate and signal.correlate switch arrays if the second array is longer.

...

...

...

signal.correlate([1, 2, 0, 0, 0], [0, 0, 1, 0, 0])

array([0, 0, 1, 2, 0, 0, 0, 0, 0])

...

...

...

signal.correlate([0, 0, 1, 0, 0],[1, 2, 0, 0, 0] )

array([0, 0, 0, 0, 0, 2, 1, 0, 0])

Well, you just happened to have very peculiar entries :) signal.correlate([-1, -2, -3], [1, 2, 3]) -> array([ -3, -8, -14, -8, -3]) signal.correlate([1, 2, 3], [-1, -2, -3]) -> array([ -3, -8, -14, -8, -3])

I found it a bit confusing, that every package needs it's own correlate function numpy, signal, ndimage, stsci.

Yes, it is. Ndimage and stsci are geared toward particular usages, signal and numpy correlate are the general functions (or maybe I am too biased toward signal processing - let's say signal.correlate is the matlab correlate). numpy.correlate only handles 1d arrays (what matlab xcorr does), whereas signal's one handle arbitrary dimension (what matlab image toolbox imfilter does). It is already quite hard to handle arbitrary dimension, but doing so in a fast way is very difficult, so this justifies at least both implementations (i.e. replacing numpy.correlate with scipy.signal one is not feasible IMHO, at least not with the current implementation). ndimage's one is much more complete, by handling many boundaries conditions (zero and constant padding, mirroring and repeating). This makes things even slower without specialized codepaths. I find it annoying that numpy.correlate and scipy.signal.correlate do not use the same defaults. Ideally, I would like that for every supported input x/y, np.correlate(x, y) == scipy.signal.correlate(x, y) (precision problems aside), but I guess this won't change. I should also note that my replacement for scipy correlate (not committed yet) is based on a new "neighborhood iterator" I have created for that purpose, so the code is much easier to follow I believe (without being much slower than the current code). This should enables implementation of faster codepaths for common cases in scipy.signal.correlate much easier (scipy.signal.correlate is often unusable for very large arrays because it is either too memory hungry or too slow - we should care about the cases where fft-based convolution is not usable, though, as fft is almost always the way to go for large arrays of comparable size). cheers, David

josef.pktd@gmail.com wrote:

On Sun, May 24, 2009 at 7:38 AM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:

...

josef.pktd@gmail.com wrote:

...

On Sun, May 24, 2009 at 6:16 AM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:

...

Hi,

I have taken a look at the correlate function in scipy.signal. There are several problems with it. First, it is wrong on several accounts: - It assumes that the correlation of complex numbers corresponds to complex multiplication, but this is not the definition followed by most textbooks, at least as far as signal processing is concerned. - More significantly, it is wrong with respect to the ordering: it assumes that correlate(a, b) == correlate(b, a), which is not true in general.

I don't see this in the results. There was recently the report on the mailing list that np.correlate and signal.correlate switch arrays if the second array is longer.

...

...

...

signal.correlate([1, 2, 0, 0, 0], [0, 0, 1, 0, 0])

array([0, 0, 1, 2, 0, 0, 0, 0, 0])

...

...

...

signal.correlate([0, 0, 1, 0, 0],[1, 2, 0, 0, 0] )

array([0, 0, 0, 0, 0, 2, 1, 0, 0])

Well, you just happened to have very peculiar entries :)

signal.correlate([-1, -2, -3], [1, 2, 3]) -> array([ -3, -8, -14, -8, -3])

signal.correlate([1, 2, 3], [-1, -2, -3]) -> array([ -3, -8, -14, -8, -3])

One of your arrays is just the negative of the other, and correlate is the same in this case. For other cases, the results differ

Gr, you're right of course :) But it still fails for arrays where the second argument has any dimension which is larger than the first one. I don't know if the assumption is relied on in the C implementation (the arrays are inverted in the C code in that case - even though they seem to be already inverted in python).

I looked at it only for examples to calculate auto-correlation and cross-correlation in time series, and had to try out to see which version works best.

Yes, it depends. I know that correlate is way too slow for my own usage in speech processing, for example for linear prediction coding. As only a few lags are necessary, direct implementation is often faster than FFT one - I have my own straightfoward autocorrelation in scikits.talkbox; I believe matlab xcorr (1d correlation) always uses the FFT. I know other people have problems with the scipy.signal correlate as well for large arrays (the C code does a copy if the inputs are not contiguous, for example - my own code using iterators should not need any copy of inputs).

Are the convolve in all cases compatible, identical (through delegation) to correlate?

I don't think so - I don't think convolution uses the conjugate for complex values. But I don't know any use of complex convolution, although I am sure there is. Correlation is always defined with the complex conjugate of the second argument AFAIK. For real cases, convolutions should always be implementable as correlation, at least when considering 0 padding for boundaries. There may be problem for huge arrays, though - doing convolution from correlation without using copies while staying fast may not always be easy, but I have never tried to do so. cheers, David

josef.pktd@gmail.com wrote:

I read somewhere in the description for one of the correlate that the code is faster if the first array is longer than the second.

Yes, but that's an implementation detail - it should not affect the result. The speed difference may come from zero padding handling. But if copies are acceptable, correlate(a, b) and correlate(b, a) can be obtained from 'inverted index', i.e., ignoring boundaries effects correlate(a, b) == correlate(b, a)[::-1, ::-1] (for rank 2 arrays). This may be much faster. FWIW, matlab is much slower at doing xcorr2(a, b) than xcorr2(b, a) if b has significantly smaller number of items.

I'm new to this, but following the discussion on the cython list, I thought that to have fast c code you always need to copy to a contiguous array.

It depends on the computation. The idea is that generally, it is really slow to do "memory jumps" on common CPU, especially when the data cannot fit the cache anymore. Doing operations one item after the other in memory is much faster than doing it 'randomly' (for a real explanation, you may consult http://people.redhat.com/drepper/cpumemory.pdf). For things like item/item computation (cos, sin, or things which are processed axis per axis), copy + contiguity is often the way to go, because this way, the items you need for computation are near in term of memory location. Now, for correlation or convolution, it is more complex, because you need to do the operation for a neighborhood. The whole problem is that the neighborhood is defined in terms of the image (spacial locality, the multi-dimensional indexes are close to each other), not in terms of memory. For rank one array, it is the same for contiguous arrays. For 2d arrays, it is already a bit difficult, but for rank 3 or rank 4, it is much more complex: for any point in a 3d image, for example, the nearest voxels can be pretty far in terms of memory. To match spatial locality to memory locality, you would need a more complex memory model than what numpy can offer out of the box, I think. In my rewriting of the correlate function, I implemented things in term of a 'neighborhood' iterator (it looks like this is nothing new - some package like itk already had this for years). The resulting code is simpler, because the complexity is hidden in the iterator, and speed is approximately the same. We are quite slower than matlab, though (a factor of 3-4 compared to imfilter). The code is there: http://github.com/cournape/scipy3/tree/clean_correlate

Are there performance differences between calling it once with a huge array versus calling it very often (eg. for sequential "online" computation) for medium sized arrays.

What is huge in your talk box applications?

In my own case, it not huge or even big. For many audio tasks, you need to compute things on a frame-per-frame basis (for example, splitting a big signal into overlapping segments of say 512 samples, and do 1d autocorrelation of this). So the usual trick is to split the big signal into a matrix where each row or column is one frame, and call a function which processes in one axis. None of the correlation function in numpy/scipy has this AFAIK. Some people sent me code for talkbox so that it can handle correlation with very big arrays, but it is one or two dimensions only. cheers, David