Correlation using FFT

Dear All,

Can someone describe the algorithm for performing correlation analysis
between two signals using fft. I tried using the simple correlation
algorithm but as my signals are quite long (upto 45 seconds of audio), as
you can imagine, the calculation takes for ever.

I know there is a quicker technique to do this by transforming the signals
into frequency spectrum using fft, but I am not sure of the exact algorithm.
I do have the advantage that both of the signals are the exact same length.

BTW, I am using kissfft as the fft library for my work.

If anyone can help me with this, or point to me a suitable link, I would be
very greatful.

Many Thanks,

RG

rg wrote:
> Dear All,
> 
> Can someone describe the algorithm for performing correlation analysis
> between two signals using fft. I tried using the simple correlation
> algorithm but as my signals are quite long (upto 45 seconds of audio), as
> you can imagine, the calculation takes for ever.
> 
> I know there is a quicker technique to do this by transforming the signals
> into frequency spectrum using fft, but I am not sure of the exact algorithm.
> I do have the advantage that both of the signals are the exact same length.
> 
> BTW, I am using kissfft as the fft library for my work.
> 
> If anyone can help me with this, or point to me a suitable link, I would be
> very greatful.
> 
> Many Thanks,
> 
> RG
> 
> 

If you are not already, you should get familiar with the concepts of 
circular convoultion and overlap-add filtering.
Many DSP texts cover this. Also, pages 6 & 7 of
http://www.borgerding.net/comp.dsp/borgerding_fastconvfilt.pdf
cover overlap-add.

The algorithm is probably more involved than I can describe quickly, but 
I'll try nonetheless.

You need to do two ffts per buffer. The buffer length is nfft minus 
twice the number of lags.  So if nlags = 128, you can use nfft = 1024 
and consume 768 samples each buffer.  One of the signals must have 
overlapping input, the other should be zero padded. I think half the 
zero padding needs to be in front, half at back; but I can't recall exactly.

Summation across buffers can be done in the frequency domain, postponing 
the inverse fft until an answer is required (i.e. only once).

A great book on the technique (as well as many others) is Richard 
Blahut's "Fast Algorithms for Digital Signal Processing".  It is out of 
print, but you can probably find it used.  Amazon.com shows 3 in the 
$65-$75 range.

-- Mark

P.S.  Glad to hear you are using kissfft. If you are willing to "give 
back", we can put your cross-correlation code into the ./tools/ section.

Hi Mark,

Thank you for your reply. It is certainly more complicated then I thought it
would be, but I will figure it out. Thanks also for the link to your
document, I will read through it.
As for the cross correlation function, I would be more then happy to provide
the code for it when I have it working. Mind you, I write my code in C++,
STL-compliant style but it should be easy to make that C compliant.

Many Thanks,

RG

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:5I2Yc.242341$fv.83510@fe2.columbus.rr.com...
> rg wrote:
> > Dear All,
> >
> > Can someone describe the algorithm for performing correlation analysis
> > between two signals using fft. I tried using the simple correlation
> > algorithm but as my signals are quite long (upto 45 seconds of audio),
as
> > you can imagine, the calculation takes for ever.
> >
> > I know there is a quicker technique to do this by transforming the
signals
> > into frequency spectrum using fft, but I am not sure of the exact
algorithm.
> > I do have the advantage that both of the signals are the exact same
length.
> >
> > BTW, I am using kissfft as the fft library for my work.
> >
> > If anyone can help me with this, or point to me a suitable link, I would
be
> > very greatful.
> >
> > Many Thanks,
> >
> > RG
> >
> >
>
> If you are not already, you should get familiar with the concepts of
> circular convoultion and overlap-add filtering.
> Many DSP texts cover this. Also, pages 6 & 7 of
> http://www.borgerding.net/comp.dsp/borgerding_fastconvfilt.pdf
> cover overlap-add.
>
> The algorithm is probably more involved than I can describe quickly, but
> I'll try nonetheless.
>
> You need to do two ffts per buffer. The buffer length is nfft minus
> twice the number of lags.  So if nlags = 128, you can use nfft = 1024
> and consume 768 samples each buffer.  One of the signals must have
> overlapping input, the other should be zero padded. I think half the
> zero padding needs to be in front, half at back; but I can't recall
exactly.
>
> Summation across buffers can be done in the frequency domain, postponing
> the inverse fft until an answer is required (i.e. only once).
>
> A great book on the technique (as well as many others) is Richard
> Blahut's "Fast Algorithms for Digital Signal Processing".  It is out of
> print, but you can probably find it used.  Amazon.com shows 3 in the
> $65-$75 range.
>
> -- Mark
>
> P.S.  Glad to hear you are using kissfft. If you are willing to "give
> back", we can put your cross-correlation code into the ./tools/ section.

"rg" <rg1117@hotmail.com> wrote in message news:<cgqf8r$hrv$1@news.freedom2surf.net>...
> Hi Mark,
> 
> Thank you for your reply. It is certainly more complicated then I thought it
> would be, but I will figure it out. Thanks also for the link to your
> document, I will read through it.

Here's a two-buffer octave/matlab script that demonstrates the principles:

nlags=8;
x=randn(32,1)+j*randn(32,1);
y=randn(32,1)+j*randn(32,1);

zpad=zeros(nlags,1);
%first buffer
X1=fft([zpad; x(1:24)] );
Y1=fft([zpad; y(1:16);zpad] );

%second buffer
X2=fft( [ x(9:16); x(17:32);zpad ] );
Y2=fft( [ zpad;    y(17:32);zpad] );

xcpw = conj( ifft( X1 .* conj( Y1 ) + X2 .* conj( Y2 ) ) );
xcp = [ xcpw(end-nlags+1:end);xcpw(1:nlags+1) ];

xc = xcorr(x,y,nlags);
snr = 10*log10( sum(abs(xc).^2)/sum(abs(xcp-xc).^2) )

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:5d6e06ef.0408281541.1577680e@posting.google.com...
> "rg" <rg1117@hotmail.com> wrote in message
news:<cgqf8r$hrv$1@news.freedom2surf.net>...
> > Hi Mark,
> >
> > Thank you for your reply. It is certainly more complicated then I
thought it
> > would be, but I will figure it out. Thanks also for the link to your
> > document, I will read through it.
>
> Here's a two-buffer octave/matlab script that demonstrates the principles:
>
> nlags=8;
> x=randn(32,1)+j*randn(32,1);
> y=randn(32,1)+j*randn(32,1);
>
> zpad=zeros(nlags,1);
> %first buffer
> X1=fft([zpad; x(1:24)] );
> Y1=fft([zpad; y(1:16);zpad] );
>
> %second buffer
> X2=fft( [ x(9:16); x(17:32);zpad ] );
> Y2=fft( [ zpad;    y(17:32);zpad] );
>
> xcpw = conj( ifft( X1 .* conj( Y1 ) + X2 .* conj( Y2 ) ) );
> xcp = [ xcpw(end-nlags+1:end);xcpw(1:nlags+1) ];
>
> xc = xcorr(x,y,nlags);
> snr = 10*log10( sum(abs(xc).^2)/sum(abs(xcp-xc).^2) )

Just a note - that if you are estimating time-delays using cross correlation
then a more refined method is necsessary where you weight the cross-spectral
density and then inverse FFT. If the weighting function is unity then you
are back to ordinary cross correlation. To get Cross PSD you can do this

Sxy(i)=beta*Sxy(i-1)+(1-beta)*X(i)*Y(i)

where X and Y are the FFTs of the two signals. Y needs to be the conjugate
in fact (or X can be conjugated instead). Beta is a forgetting factor
(0<beta<1) and i is the frame number. Once you have PSD (or
cross-periodogram as it is better known) you then do an inverse FFT. The
generalised method is better when the noises are non-white (your case of
course).Then you need a weighting factor based on coherence (Hanan-Thomson
or SCOT method - and so on in the literature)

Tom

Tom wrote:
> "Mark Borgerding" <mark@borgerding.net> wrote in message
> news:5d6e06ef.0408281541.1577680e@posting.google.com...
> 
>>"rg" <rg1117@hotmail.com> wrote in message
> 
> news:<cgqf8r$hrv$1@news.freedom2surf.net>...
> 
>>>Hi Mark,
>>>
>>>Thank you for your reply. It is certainly more complicated then I
> 
> thought it
> 
>>>would be, but I will figure it out. Thanks also for the link to your
>>>document, I will read through it.
>>
>>Here's a two-buffer octave/matlab script that demonstrates the principles:
>>
>>nlags=8;
>>x=randn(32,1)+j*randn(32,1);
>>y=randn(32,1)+j*randn(32,1);
>>
>>zpad=zeros(nlags,1);
>>%first buffer
>>X1=fft([zpad; x(1:24)] );
>>Y1=fft([zpad; y(1:16);zpad] );
>>
>>%second buffer
>>X2=fft( [ x(9:16); x(17:32);zpad ] );
>>Y2=fft( [ zpad;    y(17:32);zpad] );
>>
>>xcpw = conj( ifft( X1 .* conj( Y1 ) + X2 .* conj( Y2 ) ) );
>>xcp = [ xcpw(end-nlags+1:end);xcpw(1:nlags+1) ];
>>
>>xc = xcorr(x,y,nlags);
>>snr = 10*log10( sum(abs(xc).^2)/sum(abs(xcp-xc).^2) )
> 
> 
> Just a note - that if you are estimating time-delays using cross correlation
> then a more refined method is necsessary where you weight the cross-spectral
> density and then inverse FFT. If the weighting function is unity then you
> are back to ordinary cross correlation. To get Cross PSD you can do this
> 
> Sxy(i)=beta*Sxy(i-1)+(1-beta)*X(i)*Y(i)

Something doesn't seem right.
I can't see any value of beta that would make the above equivalent to 
normal (i.e. nonleaky) cross-correlation.
i.e.
Sxy(i) = Sxy(i-1)  +  X(i)*Y(i)

To make a "leaky integrator" fast cross correlator,  I would  omit the 
second coefficient in your formula, 1-beta.
Leaving,
Sxy(i) = beta*Sxy(i-1)  +  X(i)*Y(i)

This makes the boundary cases nice and clean:
beta = 0 -- forget everything , use current block only
beta = 1 -- forget nothing , use entire signal


-- Mark

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:4135c941$1@news.xetron.com...
> To make a "leaky integrator" fast cross correlator,  I would  omit the
> second coefficient in your formula, 1-beta.
> Leaving,
> Sxy(i) = beta*Sxy(i-1)  +  X(i)*Y(i)
>
> This makes the boundary cases nice and clean:
> beta = 0 -- forget everything , use current block only
> beta = 1 -- forget nothing , use entire signal
>
>
> -- Mark
>
If you do what you suggest you get a dc gain (ie when z=1) which is
1/(1-beta) rather than unity.
Take z-transforms and the TF is

(1-beta)/(1-betaz^-1) ....

I see your argument - you are trying to get back to a pure integrator when
beta=1 but pure integrators are not a good idea in open-loop - any slight
dc-offset and off they go for a walk.Best results are obtained with beta=05
up to 0.9 (ish!).

Tom

Tom wrote:
> "Mark Borgerding" <mark@borgerding.net> wrote in message
> news:4135c941$1@news.xetron.com...
> 
>>To make a "leaky integrator" fast cross correlator,  I would  omit the
>>second coefficient in your formula, 1-beta.
>>Leaving,
>>Sxy(i) = beta*Sxy(i-1)  +  X(i)*Y(i)
>>
>>This makes the boundary cases nice and clean:
>>beta = 0 -- forget everything , use current block only
>>beta = 1 -- forget nothing , use entire signal
>>
>>
>>-- Mark
>>
> 
> If you do what you suggest you get a dc gain (ie when z=1) which is
> 1/(1-beta) rather than unity.
> Take z-transforms and the TF is
> 
> (1-beta)/(1-betaz^-1) ....
> 
> I see your argument - you are trying to get back to a pure integrator when
> beta=1 but pure integrators are not a good idea in open-loop - any slight
> dc-offset and off they go for a walk.Best results are obtained with beta=05
> up to 0.9 (ish!).
> 
> Tom

I agree pure integration is a bad idea for endless input, but that was 
not the problem put forth.
The OP didn't mention anything about an open loop.

FWIW, I don't think the term "unity gain" is very meaningful when 
applied to a single buffer answer in cross-correlation.

Let's hop a little further down this bunny trail ...

Both systems' transfer functions contain a single pole on the real axis 
with magnitude beta.

I suggest that the unity gain created by the (1-beta) term causes more 
problems than
it solves.  It certainly declaws the unstable pole when beta == 1 by 
setting the gain to zero.
Unfortunately, that happens to be is the useful case of 
cross-correlation of two complete sequences.

To have the best of both worlds, I'd split beta and 1-beta up into two 
gains: beta and alpha.
	y(i) = beta*y(i-1)  +  alpha*x(i)
For the case when unity gain is desired,  alpha = 1-beta.  If 
integration is the goal, then alpha = beta = 1

Perhaps it comes down to personal preference.
I prefer one algorithm that does two things even if it is slightly more 
complicated, rather than needing two algorithms.

In any case, a gain is usually easy to slip in someplace computationally 
convenient.

-- Mark

"Mark Borgerding" <mark@borgerding.net> wrote in message
news:41374d72$1@news.xetron.com...
> Tom wrote:
> > "Mark Borgerding" <mark@borgerding.net> wrote in message
> > news:4135c941$1@news.xetron.com...
> >
> >>To make a "leaky integrator" fast cross correlator,  I would  omit the
> >>second coefficient in your formula, 1-beta.
> >>Leaving,
> >>Sxy(i) = beta*Sxy(i-1)  +  X(i)*Y(i)
> >>
> >>This makes the boundary cases nice and clean:
> >>beta = 0 -- forget everything , use current block only
> >>beta = 1 -- forget nothing , use entire signal
> >>
> >>
> >>-- Mark
> >>
> >
> > If you do what you suggest you get a dc gain (ie when z=1) which is
> > 1/(1-beta) rather than unity.
> > Take z-transforms and the TF is
> >
> > (1-beta)/(1-betaz^-1) ....
> >
> > I see your argument - you are trying to get back to a pure integrator
when
> > beta=1 but pure integrators are not a good idea in open-loop - any
slight
> > dc-offset and off they go for a walk.Best results are obtained with
beta=05
> > up to 0.9 (ish!).
> >
> > Tom
>
> I agree pure integration is a bad idea for endless input, but that was
> not the problem put forth.
> The OP didn't mention anything about an open loop.
>
> FWIW, I don't think the term "unity gain" is very meaningful when
> applied to a single buffer answer in cross-correlation.
>
> Let's hop a little further down this bunny trail ...
>
> Both systems' transfer functions contain a single pole on the real axis
> with magnitude beta.
>
> I suggest that the unity gain created by the (1-beta) term causes more
> problems than
> it solves.  It certainly declaws the unstable pole when beta == 1 by
> setting the gain to zero.
> Unfortunately, that happens to be is the useful case of
> cross-correlation of two complete sequences.
>
> To have the best of both worlds, I'd split beta and 1-beta up into two
> gains: beta and alpha.
> y(i) = beta*y(i-1)  +  alpha*x(i)
> For the case when unity gain is desired,  alpha = 1-beta.  If
> integration is the goal, then alpha = beta = 1
>
> Perhaps it comes down to personal preference.
> I prefer one algorithm that does two things even if it is slightly more
> complicated, rather than needing two algorithms.
>
> In any case, a gain is usually easy to slip in someplace computationally
> convenient.
>
> -- Mark
>
The idea comes from exponential smoothing in time-series analysis
http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc431.htm

If you don't put in the (1-beta) part then you get an offset ie you have to
re-scale afterwards to get the right
answer.
You will find that integrators (pure) are rarely if ever used on open loop
(as we have here).
For example the LMS algorithm has integrators in it - as does the Kalman
Filter and many such algorithms

w(k+1)=w(k)+2*mu*X(k)*e(k)

but it has an error term ie feedback to keep it in line!It does not matter
whether integrators are analogue or digital, in open loop they are
troublesome.Try putting beta=1 and see if it works.Or simpler still try
integrating a pure sine-wave.Slightest dc and we are in trouble.

Tom

Tom wrote:
> If you don't put in the (1-beta) part then you get an offset ie you have to
> re-scale afterwards to get the right
> answer.

On the contrary, it is impossible to get the right answer if you DO put 
in the 1-beta part.

... if the question being answered is the one the original poster asked: 
how to correlate two sequences using the fft.

Or perhaps I am also guilty of assuming too much.   Perhaps the OP 
wanted the single-valued correlation : cov(x,y) / sqrt( var(x)*var(y) )
But I don't think so.  Covariance and variance operations are already 
linear complexity.  I can't see how FFTs would make them any faster. 
I'm pretty sure the OP wanted cross-correlation when he asked for 
correlation.

There is a very specific definition for cross-correlation.
	See http://mathworld.wolfram.com/Cross-Correlation.html
	or http://cnx.rice.edu/content/m10686/latest/
Notice the pure integration and infinite bounds.

The formula you posted may have uses for continuous (open loop) 
processing, but it is NOT cross-correlation.

I can only assume the question for your "right answer" is "how do I get 
a time-decaying approximation of cross-correlation".  That question was 
never asked.  You made a good suggestion that was certainly topical, 
considering more people read a thread than just the few persons writing 
it. But the OP had signals "up to 45 seconds long", far from infinite.

-- Mark