[PD] Cross correlation using FFT

[Charles Henry]

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 at postlude.co.uk> wrote:
> On Thu, 6 Apr 2006 18:54:55 +0000
> Jamie Bullock <jamie at 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 at iem.at mailing list
> UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
>