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/b...stconvfilt.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" wrote in message
...
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/b...stconvfilt.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" wrote in message ...
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" wrote in message
om...
"rg" wrote in message
...
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
(0beta1) 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" wrote in message
om...

"rg" wrote in message

...

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" wrote in message
...
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" wrote in message
...

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" wrote in message
...
Tom wrote:
"Mark Borgerding" wrote in message
...

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/handb...on4/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

"Mark Borgerding" wrote in message
...
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/know...anshuang-1.pdf

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

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