# Correlation using FFT

Last response: in Home Audio

rg

August 28, 2004 6:52:47 PM

Archived from groups: comp.dsp,rec.audio.tech,sci.math.num-analysis (More info?)

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

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

More about : correlation fft

Anonymous

August 28, 2004 8:49:05 PM

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_fastconvf...

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

August 28, 2004 10:25:29 PM

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_fastconvf...

> 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.

Anonymous

August 28, 2004 10:25:30 PM

"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) )

Tom

September 1, 2004 8:12:42 PM

"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

Anonymous

September 1, 2004 8:12:43 PM

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

Tom

September 3, 2004 12:00:58 AM

"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

Anonymous

September 3, 2004 12:00:59 AM

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

Tom

September 3, 2004 8:50:36 PM

"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/pm...

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

Anonymous

September 3, 2004 8:50:37 PM

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

Tom

September 5, 2004 12:31:02 AM

"Mark Borgerding" <markab@xetron.com> wrote in message

news:413871c0$1@news.xetron.com...

> 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

> ]

It is certainly Cross Correlation. In fact the method is more general (with

modification) - I mentioned the Hannan-Thomson algorithm and the SCOT method

for instance.

In its simpler form it is just the Wiener-Kinchen theorem ie the Fourier

Transform of cross-corr is power spectral density.

The basic idea is found here

http://www.paper.edu.cn/scholar/known/qiutianshuang/qiu...

though it is much older.

You are confusing the pure integration in smoothing the cross-periodograms

with the pure integration of getting a Cross-Corr (as you point out). The

method I quoted is only to smooth the cross periodogram (PSD estimate).You

don't need pure integration as this happens when you do an inverse FFT ie

the Wiener Kinchen theorem again.

In fact you don't need to smooth the cross-periodograms at all but you would

get quite a noisy estimate.

Also ordinary cross-corr whether you do it in the freq domain or the

ordinary method has drawbacks with certain problems - particulalry

estimating time-delays (or multple time-delays) and this is why generalised

cross-corr is necessary.

Hope this helps you to understand.

Tom

Anonymous

September 7, 2004 11:41:46 AM

Mark Borgerding wrote:

....

> 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

I'm one other reading the thread...

maybe my contribution is far off the original goal of the OP.

Just to satisfy my interest (and enhance my knowledge...), I would

like to put in two more general questions:

1) given two such signals of 45 seconds, which are equal.

now cut 5s from the beginning of the first and paste it to the

end of it .

then correlating this modified signal with the other unmodified.

I reckon that the frequency content is the same except on the

borders where the signal has been cut/pasted.

Therefore I would expect similar FFT results and probably a high

degree of correlation.

My question: does correlation detect the shift of most of the

signal and show good correlation?

Or would the result be a low correlation because signal contents

don't match at any point?

2) Given, I want to compare a signal (endless stream) with a certain

pattern (short piece of signal).

I expect that the signal contains pieces which are similar to the

pattern, but differ either

* in amplitude or

* in width or

* in both.

Which method would be chosen to find the occurrences?

Is cross-correlation the right approach?

Or will it fail in one of the cases?

Remember: it's not a current task which bothers me,

it's only that I'm curious

(maybe I'll need the answer soon, nevertheless ...)

so it's enough to give me a direction.

Bernhard

Read discussions in other Home Audio categories

!