Hi,
Has anyone implemented cross correlation in PD using rifft~? I am confident with performing cross correlation in the time domain, but am a little confused about how the equivalent operation can be conducted using FFTs. So far I understand the algorithm to be:
- Zero pad the signal vectors to twice the vector size
- Take their ffts: fft(x) fft(y)
- Take the complex conjugate of fft(x) and multiply this by fft(y)
- Take the inverse fft of the result.
Whatever I try doesn't seem to give the correct output. In particular - I am unsure of step 3 above. Can anyone give some clarification or an example patch.
TIA,
Jamie
On Mon, 3 Apr 2006, Jamie Bullock wrote:
I'll summarize what we said on the #dataflow chatline:
- Zero pad the signal vectors to twice the vector size
This is to prevent the result from wrapping around.
- Take their ffts: fft(x) fft(y)
- Take the complex conjugate of fft(x) and multiply this by fft(y)
Whatever I try doesn't seem to give the correct output. In particular - I am unsure of step 3 above. Can anyone give some clarification or an example patch.
there are several ways to do steps 2 and 3, because there are several identities about fft. especially, fft(fft(signal)) = reverse(signal) * blocksize. Note that reverse() here keeps sample t=0 at the same place and only reverses the rest of the block.
Then fft(fft(fft(fft(signal)))) = signal * blocksize^2
Then ifft is like fft of the reversal, and also like the conj of the fft.
Pd now can do FFT on pictures, using GridFlow's [#fft]. I mention it in case it helps you visualizing what's going on. You will need to install FFTW.
_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju | Freelance Digital Arts Engineer, Montréal QC Canada
Hi,
On Thu, 6 Apr 2006 13:35:04 -0400 (EDT) Mathieu Bouchard matju@artengine.ca wrote:
On Mon, 3 Apr 2006, Jamie Bullock wrote:
I'll summarize what we said on the #dataflow chatline:
- Zero pad the signal vectors to twice the vector size
This is to prevent the result from wrapping around.
- Take their ffts: fft(x) fft(y)
- Take the complex conjugate of fft(x) and multiply this by fft(y)
Whatever I try doesn't seem to give the correct output. In particular - I am unsure of step 3 above. Can anyone give some clarification or an example patch.
there are several ways to do steps 2 and 3, because there are several identities about fft. especially, fft(fft(signal)) = reverse(signal) * blocksize. Note that reverse() here keeps sample t=0 at the same place and only reverses the rest of the block.
Then fft(fft(fft(fft(signal)))) = signal * blocksize^2
Then ifft is like fft of the reversal, and also like the conj of the fft.
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...
Pd now can do FFT on pictures, using GridFlow's [#fft]. I mention it in case it helps you visualizing what's going on. You will need to install FFTW.
Thanks. I'll check that out.
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
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
Hi,
On Fri, 7 Apr 2006 12:38:53 -0500 "Charles Henry" czhenry@gmail.com wrote:
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
And if you weren't time reversing y, it would be:
(xy + x'y') + (xy' + x'y)i
.. as Matju suggests?
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) )
I'm afraid I don't understand your notation here. Are we still dealing with FFTs? ...and is i the subscript for the input vector array y?
Something like:
for(i = 0; i < ??? ; i++){
for(j = 1; j < N; j++)
correlation_coefficient[i] = x[j] * y[j - i];
}
.. or are x and y fft(x) and fft(y), making the above code somewhat more complex ?
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.
Did you devise that approach yourself, or is there a paper you can point me towards?
There are some problems with this, to be sure. It all depends on what you want to use it for?
I would like to try a few things like detecting proximity between two microphones, and detection of acoustic spill. Maybe also fundamental estimation by auto correlation. It just seems a useful tool.
Chuck
P.S. I don't have the external available at the moment....but I can send it later, if you'd like.
That would be really good. I think I probably understand code better than I understand maths ;-)
We can do a swap - you can find my time-domain circular cross correlation external as part of my flib library in the PD CVS under postlude. It's just a brute force approach, not especially pretty or efficient. An optional argument gives the maximum delay (in samples) between the two inputs for each correlation coefficient. Conceptually the summations are performed for -maxdelay < 0 < maxdelay, but in reality the values are 'wrapped' so that the index to the delayed signal stays within the array bounds.
best, Jamie
On Fri, 7 Apr 2006, Jamie Bullock wrote:
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 ?
no, the case when imaginaries cancel happens only when you multiply a number by its own conjugate, and then it gives the square of the magnitude:
(x + ix')(y - iy') = xy + xiy' + ix'y - ix'iy' = (xy + x'y') + (xy' + x'y)i
so in the case when x=y and x'=y',
(x + ix')(x - ix') = (xx + x'x') + (xx' + x'x)i = (x^2 + x'^2) + 0i
_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju | Freelance Digital Arts Engineer, Montréal QC Canada
On Fri, 7 Apr 2006, Mathieu Bouchard wrote:
no, the case when imaginaries cancel happens only when you multiply a number by its own conjugate, and then it gives the square of the magnitude: (x + ix')(y - iy') = xy + xiy' + ix'y - ix'iy' = (xy + x'y') + (xy' + x'y)i
That and the ordinary complex product should be abstractions, but Miller doesn't seem to use abstractions much in his patches (see his FFT tutorials...), so the two kinds of complex products are copypasted all over the place instead. Four [*~], one [+~], one [-~], and eight more wires than using an abstr. Brilliant.
_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju | Freelance Digital Arts Engineer, Montréal QC Canada
On Fri, 7 Apr 2006, Mathieu Bouchard wrote:
On Fri, 7 Apr 2006, Mathieu Bouchard wrote:
no, the case when imaginaries cancel happens only when you multiply a number by its own conjugate, and then it gives the square of the magnitude: (x + ix')(y - iy') = xy + xiy' + ix'y - ix'iy' = (xy + x'y') + (xy' + x'y)i
That and the ordinary complex product should be abstractions, but Miller doesn't seem to use abstractions much in his patches (see his FFT tutorials...), so the two kinds of complex products are copypasted all over the place instead. Four [*~], one [+~], one [-~], and eight more wires than using an abstr. Brilliant.
D'OH, typo:
(x + ix')(y - iy') = xy - xiy' + ix'y - ix'iy' = (xy + x'y') + (x'y - xy')i
_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju | Freelance Digital Arts Engineer, Montréal QC Canada
Sorry it took me so long to send the external and code..... It just got me thinking that I was doing something really wrong here, but I think I've got it straightened out now.
The implementation here does a pretty good job of using xcorr~ in overlapping blocks, like block 256 2. So, we have the two signals, as described before sig1 |////|////| 2N non-zero samples sig2 |----//|//----| N non-zero samples, in the middle We can do a running calculation of the cross-correlation, using overlapping blocks. (Don't reblock the output, it won't work that way.) |////|////| |----//|//----| |--------|////|////| |--------|----//|//----|
And the result is that when we add up the results of xcorr across blocks we get a running calculation of the cross correlation. So, I also created two others to help here. intbb~ which integrates (stably) between blocks (which was no doubt unnecessary) and maxm~ which picks out the maximum value, and its location in the vector between bounds, e.g. [maxm~ -256 256] picks out the location of the maximum value, closest to the "origin" of the vector, in either the last 256 samples of xcorr~ or the first 256 samples of xcorr~
again sorry for the long delay chuck
On 4/7/06, Mathieu Bouchard matju@artengine.ca wrote:
On Fri, 7 Apr 2006, Mathieu Bouchard wrote:
On Fri, 7 Apr 2006, Mathieu Bouchard wrote:
no, the case when imaginaries cancel happens only when you multiply a number by its own conjugate, and then it gives the square of the magnitude: (x + ix')(y - iy') = xy + xiy' + ix'y - ix'iy' = (xy + x'y') + (xy' + x'y)i
That and the ordinary complex product should be abstractions, but Miller doesn't seem to use abstractions much in his patches (see his FFT tutorials...), so the two kinds of complex products are copypasted all over the place instead. Four [*~], one [+~], one [-~], and eight more wires than using an abstr. Brilliant.
D'OH, typo:
(x + ix')(y - iy') = xy - xiy' + ix'y - ix'iy' = (xy + x'y') + (x'y - xy')i
_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju | Freelance Digital Arts Engineer, Montréal QC Canada
PD-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list
-- Charles Zachary Henry
anti.dazed.med Med student who needs a Mickey's
-
Charles Henry -
Jamie Bullock -
Mathieu Bouchard