Forums

Cross correlation via Hartley Transform

Started by pete26 September 8, 2010
Hi ... 

I am triing to get a cross correlation between two or (later 16) audio
channels with the hartley transform (I am doing that in Labview with FFT
comparison). But the results don't fit:

-> my source: http://rsbweb.nih.gov/nih-image/download/n...docs/thesis.pdf
on page 12. 

So I am doing the Hartley transform on both signals H1(f) and H2(f) ... to
get the cross correlation -> iFHT(1/2[ H1(f) * H2(f) + H1(-f) * H2(-f) +
H1(f) * H2(-f) - H1(-f) * H2(f) ]). 

I am really appreciated about every tip why this doesnt fit with the result
per FFT ? 

Thank you very much 
Peter


I am sorry the topic is twice in the forum ... can anyone of the moderators
delete this one.

Thank you very much !


>Hi ... > >I am triing to get a cross correlation between two or (later 16) audio >channels with the hartley transform (I am doing that in Labview with FFT >comparison). But the results don't fit: > >-> my source:
http://rsbweb.nih.gov/nih-image/download/n...docs/thesis.pdf
>on page 12. > >So I am doing the Hartley transform on both signals H1(f) and H2(f) ...
to
>get the cross correlation -> iFHT(1/2[ H1(f) * H2(f) + H1(-f) * H2(-f) + >H1(f) * H2(-f) - H1(-f) * H2(f) ]). > >I am really appreciated about every tip why this doesnt fit with the
result
>per FFT ? > >Thank you very much >Peter > > >
On Sep 8, 8:31&#2013266080;am, "pete26" <p2222@n_o_s_p_a_m.sbox.tugraz.at> wrote:
> I am sorry the topic is twice in the forum ... can anyone of the moderators > delete this one. > > Thank you very much ! > > >Hi ... > > >I am triing to get a cross correlation between two or (later 16) audio > >channels with the hartley transform (I am doing that in Labview with FFT > >comparison). But the results don't fit: > > >-> my source: > > http://rsbweb.nih.gov/nih-image/download/n...docs/thesis.pdf > > > > >on page 12. > > >So I am doing the Hartley transform on both signals H1(f) and H2(f) ... > to > >get the cross correlation -> iFHT(1/2[ H1(f) * H2(f) + H1(-f) * H2(-f) + > >H1(f) * H2(-f) - H1(-f) * H2(f) ]). > > >I am really appreciated about every tip why this doesnt fit with the > result > >per FFT ? > > >Thank you very much > >Peter- Hide quoted text - > > - Show quoted text -
You need to fix your link (didn't work for me). It's difficult to comment without it. Are you doing circular cross-correlation or zero padding to do linear? (From your other post): How are you 'flipping' those (-f) terms? Kevin McGee
>You need to fix your link (didn't work for me). It's difficult to >comment without it.
Hi, I tried to fix the link but it doesnt work on this way -> can be found by "optimized fast hartley transform for the mc68000" on page 12. I am just doing the hartley transform on the signals and then with flipping (building up the reverse vector,like "fliplr" in matlab) the vector i get H(-f). Because of using circular or linear correlation -> I am using circular cross correlation (in fact I am not really sure about it) and I am just triing to get the same result as with the FFT by doing muliplikation of the spectra ifft(fft(sig1).*fft(sig2)). Maybe I should try the formular for circular correlation on page 14 but I am not sure how to build up H(N-f) is it just the flipped version and shifted to the right with N ? Thank you for your effort Peter
> >Are you doing circular cross-correlation or zero padding to do >linear? (From your other post): How are you 'flipping' those (-f) >terms? > >Kevin McGee >
On Sep 9, 5:56&#2013266080;am, "pete26" <p2222@n_o_s_p_a_m.sbox.tugraz.at> wrote:
> >You need to fix your link (didn't work for me). &#2013266080;It's difficult to > >comment without it. > > Hi, I tried to fix the link but it doesnt work on this way -> can be found > by "optimized fast hartley transform for the mc68000" on page 12. > > I am just doing the hartley transform on the signals and then with flipping > (building up the reverse vector,like "fliplr" in matlab) the vector i get > H(-f). > > Because of using circular or linear correlation -> I am using circular > cross correlation (in fact I am not really sure about it) and I am just > triing to get the same result as with the FFT by doing muliplikation of the > spectra ifft(fft(sig1).*fft(sig2)). Maybe I should try the formular for > circular correlation on page 14 but I am not sure how to build up H(N-f) is > it just the flipped version and shifted to the right with N ? > > Thank you for your effort > Peter > > > > > > >Are you doing circular cross-correlation or zero padding to do > >linear? &#2013266080;(From your other post): How are you 'flipping' those (-f) > >terms? > > >Kevin McGee- Hide quoted text - > > - Show quoted text -
Ok, I've got that reference. But I don't use Matlab. When doing circular cross-correlations in C or C++, I get the FFT of both signals (N points each), then conjugate one of them, multiply FFT results, point by point, and inverse transform. The cross-correlation result gives me the circular cross-correlation of the two signals. If there is a peak in the result, it tells me the (integer) time delay between the signals. You also often have to interpolate the results if you want a more precise (non-integer) delay. As to which result to conjugate - it depends on how I want to interpret a +/- delay. If I want to do a linear cross-correlation, I zero pad the data such that I have 2N points for each in the time domain. That way I don't get any 'wrap around effect' (see link below) http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/correlate/ The H[N-k] (or H(-f) if you prefer) should be 'flipped' the following way: take all the data from points 1 to N-1 and reverse them - leave the zero point alone (how do you 'flip' an x(0) point? - simple - you don't). So, for example the data sequence 2, 7, 9, 4, 5, 1, 6, 8 would become: 2, 8, 6, 1, 5, 4, 9, 7. I think the formula on p. 14 is correct, excepting the obvious typo of H1(-k), which should be H1(N-k). Make sure you're getting the FFT version right first. If you want, you might even try the standard formula for cross-correlation as shown in the above link (now you've got three ways of computing things). Kevin McGee
>On Sep 9, 5:56=A0am, "pete26" <p2222@n_o_s_p_a_m.sbox.tugraz.at> wrote: >> >You need to fix your link (didn't work for me). =A0It's difficult to >> >comment without it. >> >> Hi, I tried to fix the link but it doesnt work on this way -> can be
foun=
>d >> by "optimized fast hartley transform for the mc68000" on page 12. >> >> I am just doing the hartley transform on the signals and then with
flippi=
>ng >> (building up the reverse vector,like "fliplr" in matlab) the vector i
get
>> H(-f). >> >> Because of using circular or linear correlation -> I am using circular >> cross correlation (in fact I am not really sure about it) and I am just >> triing to get the same result as with the FFT by doing muliplikation of
t=
>he >> spectra ifft(fft(sig1).*fft(sig2)). Maybe I should try the formular for >> circular correlation on page 14 but I am not sure how to build up H(N-f)
=
>is >> it just the flipped version and shifted to the right with N ? >> >> Thank you for your effort >> Peter >> >> >> >> >> >> >Are you doing circular cross-correlation or zero padding to do >> >linear? =A0(From your other post): How are you 'flipping' those (-f) >> >terms? >> >> >Kevin McGee- Hide quoted text - >> >> - Show quoted text - > > >Ok, I've got that reference. But I don't use Matlab. When doing >circular cross-correlations in C or C++, I get the FFT of both signals >(N points each), then conjugate one of them, multiply FFT results, >point by point, and inverse transform. The cross-correlation result >gives me the circular cross-correlation of the two signals. If there >is a peak in the result, it tells me the (integer) time delay between >the signals. You also often have to interpolate the results if you >want a more precise (non-integer) delay. > >As to which result to conjugate - it depends on how I want to >interpret a +/- delay. > >If I want to do a linear cross-correlation, I zero pad the data such >that I have 2N points for each in the time domain. That way I don't >get any 'wrap around effect' (see link below) > >http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/correlate/ > >The H[N-k] (or H(-f) if you prefer) should be 'flipped' the following >way: take all the data from points 1 to N-1 and reverse them - leave >the zero point alone (how do you 'flip' an x(0) point? - simple - you >don't). So, for example the data sequence 2, 7, 9, 4, 5, 1, 6, 8 >would become: 2, 8, 6, 1, 5, 4, 9, 7. > >I think the formula on p. 14 is correct, excepting the obvious typo of >H1(-k), which should be H1(N-k). > >Make sure you're getting the FFT version right first. If you want, >you might even try the standard formula for cross-correlation as shown >in the above link (now you've got three ways of computing things). > >Kevin McGee >
Thank you very much for your effort ... I'll give a comment if it worked later. Peter
On Sep 8, 8:15&#2013266080;am, "pete26" <p2222@n_o_s_p_a_m.sbox.tugraz.at> wrote:
> So I am doing the Hartleytransformon both signals H1(f) and H2(f) ... to > get the cross correlation -> iFHT(1/2[ H1(f) * H2(f) + H1(-f) * H2(-f) + > H1(f) * H2(-f) - H1(-f) * H2(f) ]).
Presumably you are computing a discrete Hartley transform, not a Hartley transform. Maybe you are confused about the indexing of the discrete transform, and in particular which indices correspond to "positive" and "negative" frequencies? (This is a common source of confusion with all of the discrete Fourier-like transforms, especially since aliasing introduces some ambiguity.) The correct formula for circular convolutions with discrete Hartley transforms, equivalent to an FFT-based convolution, is given on the Wikipedia page: http://en.wikipedia.org/wiki/Discrete_Hartley_transform (Note that, contrary to common claims, there is no known computational advantage to using a fast Hartley transform [FHT] for real convolution operations compared to using an FFT algorithm specialized for real data. The latter actually requires slightly less arithmetic in theory, and in practice it is much easier to find a highly optimized real-data FFT than it is to find a highly optimized FHT, and performance of these kinds of algorithms is nowadays mostly determined by implementation issues [see e.g. http://cnx.org/content/m16336/ ].) Regards, Steven G. Johnson
Thank you for your help ... 

I could manage to get the same results as with the fft methode but I got no
real advantage of it -> in literature it is written that you will need with
hartley transform less memory because of using real values. 

But in fact its like you mentioned ... in Labview, Matlab an FFT is that
fast and optimized and  thus not to beat with the FHT in memory and speed
sense.

By the way ... I tried to use the goertzel algorithms as well but its the
same just from 10000 samples the speed of goertzel (with 10 frequency bins)
is higher than with the fft.

Soo i guess there is no faster way to calculate a cross correlation between
two signals (length is 1920 Samples) as with fft's.

Thank you very much !

Peter