Hi, I encountered this problem before when I wanted to do a 'windowed' cross-correlation function in pd too. Your algorithm is simpler than mine, I wrote an external to do it on its own.
(x + x'i)(y-y'i) = (xy + x'y') + (x'y - xy')*i
Now then, there's the problem of finding what the output means: the cross-correlation (in this example) is being done as a convolution between signal x, and time-reversed signal y.
So, the summation notation can shed some light on it:
1st: simple convolution: x * y (i) = sum( j=1, N-1: x(j) y(i-j) ) In the case of ffts , we do circular convolution, so if it's y(N) it means y(0), and y(-5)=y(N-5)
The real thing: x * rev(y) (i) = sum ( j=1, N : x(j) rev(y) (i-j) ) = sum ( j=1, N : x(j) y (j-i) )
So, the output should reflect that for each positive i , we get the signal y, shifted i samples forward in time, in product with x. For negative i (the samples N-i), we get the signal y delayed by i samples in product with x.
In my implementation, I kept the middle half of y the same, and let x be completely the same, like this:
x:|//////////| y:|-----/////-----|
Then, the cross correlation has a maximum value in either (1, N/4) where y lags behind x, or in (3*N/4, N-1) where y leads x. Then, the samples from (N/4, 3*N/4) were just garbage, left over from the fft.
There are some problems with this, to be sure. It all depends on what you want to use it for?
Chuck
P.S. I don't have the external available at the moment....but I can send it later, if you'd like.
On 4/7/06, Jamie Bullock jamie@postlude.co.uk wrote:
On Thu, 6 Apr 2006 18:54:55 +0000 Jamie Bullock jamie@postlude.co.uk wrote:
<snip> > > Please find attached, my attempt at making this into a PD patch. The array 'proper' contains a snapshot of the output of from my PD extern, which performs time domain cross-correlation, and what I believe the output should look like. > > The output to the 'output' array is clearly wrong, but I'm not sure why. > > The complex product subpatch is based on one of Miller's examples, I can't remember which. Not sure if the -1 multiply for complex conj (added by me) is implemented properly... >
I'm still trying to work this out and getting into more knots. I thought I'd strip it back to the basic maths.
Am I right in saying that the basic operation FFT(x) * conj(FFT(y)) works like this:
(x + ix')(y - iy') = (xy + x'y') + 0i
?
Jamie
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list