Reply by Mark Borgerding●September 8, 20042004-09-08
Bernhard Holzmayer wrote:
> Mark Borgerding wrote:
> I'm one other reading the thread...
> maybe my contribution is far off the original goal of the OP.
You raise some good points. Ones that should've been addressed early on.
In hand-waving terms, cross-correlation can be applied to
a) two short signals,
b) one long + one short, or
c) two long signals
For the first two scenarios, it is feasible to compute the
cross-correlation in its entirety. The last scenario is basically
defined as one where it is impractical to calculate the entire CC
sequence. The latter case is what I've been talking about.
For a), you can transform both signals into the frequency domain by zero
padded ffts (length >= N+M-1), then conjugate multiply their spectral
components to affect cross-correlation in the time domain.
Note cross-correlation is is just time-reverse, conjugate convolution.
Also note multiplication in the freq domain <==> conv in time D.
Also^2 note conjugate in the freq D is time reversal and conjugation in
time D.
To apply this idea back to the original poster: To cross-correlate two
audio signals of 45s at 44kS/s, ( roughly 2e6 samples each), an FFT of
4e6 points would be needed. Probably feasible for the OP's purpose. I
probably should've suggested this straightaway due to its simplicity.
tic;
nf = ng = 2e6;
nfft=2^nextpow2(nf+ng-1);
f=randn(nf,1);
g=randn(ng,1);
F=fft(f,nfft);
G=fft(g,nfft);
cc=ifft( F.*conj(G) ); % time-shifted cross-correlation
toc
It took less than 30s on my computer. Sorry rg.
Anyway ...
For b), you can extend the above solution to work in a buffer-by-buffer
case by applying overlap-add or overlap-save filtering, with the shorter
sequence reversed and conjugated.
For case c), it is impractical to compute the entire non-zero
cross-correlation sequence. It is easy to calculate only a range of
delay(aka lag) values.
This is where the matlab script I posted earlier comes into play.
It performs blockwise computation of only those cross-correlation values
of interest. It postpones the ifft until an answer is needed by
exploting the identity f + g <==> F + G.
> 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.
Yes. In fact for this case, they two would vary only in phase, |F1| == |F2|
There would be a linear phase shift in the frequency domain
corresponding to the constant time shift in the time domain.
(Think about linear phase(AKA constant group delay) FIR filters)
> 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?
The cross-correlation sequence would be large and positive at the lag
where t = 5s
The CC seq is maximized when the two sequences have matching spectral
phase components that match (linear scaling of either only affects the
magnitude, of the CC).
>
> 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?
This is close to case b) above.
Caveat on "width": Cross-correlation will not "stretch" a pattern to
find a match. If a feature is "stretched", while keeping its shape, its
frequency spectrum is shifted.
Perhaps there is a cepstral technique that could accomplish this task.
>
> 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
I hope I helped.
-- Mark
Reply by Bernhard Holzmayer●September 8, 20042004-09-08
Mark Borgerding wrote:
>> 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?
>
> This is close to case b) above.
>
> Caveat on "width": Cross-correlation will not "stretch" a pattern
> to
> find a match. If a feature is "stretched", while keeping its
> shape, its frequency spectrum is shifted.
>
> Perhaps there is a cepstral technique that could accomplish this
> task.
>
Thanks for taking the challenge.
Measuring the correlation with such a normalized pattern will
probably be a future task, so I'm really interested in any ideas.
What I bear in mind to solve this task has grown from the wavelet
theory, and your response reminds me of that:
If I generate a series of related patterns which differ only in
width (time), sort of correlation with every pattern should bring
up a vector which gives me a "best fit" result.
Instead of working with generic sine and/or cosine signals which
would certainly lead to something like FFT, it seems as if this
were more a task for wavelet transform, which uses short patterns.
Currently, I have no idea which would be the best approach, not
even, if the results would be acceptable for a real application,
just looking around ...
Bernhard
Reply by Bernhard Holzmayer●September 7, 20042004-09-07
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
Reply by Tom●September 4, 20042004-09-04
"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/qiutianshuang-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
Reply by Mark Borgerding●September 3, 20042004-09-03
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
Reply by Tom●September 3, 20042004-09-03
"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
Reply by Mark Borgerding●September 2, 20042004-09-02
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
Reply by Tom●September 2, 20042004-09-02
"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
Reply by Mark Borgerding●September 1, 20042004-09-01
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
Reply by Tom●September 1, 20042004-09-01
"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