Hartley Transform

Entry posted by EricBall November 5, 2009

504 views

If you do any kind of signal processing, you've probably heard of an "FFT" or Fast Fourier Transform. An FFT is an algorithm (and there are several) which calculates a Discrete Fourier Transform in less operations, typically O(N*log2(N)), where N is the number of samples. The Fourier Transform changes a time based function into a frequency based function. (The reverse is also possible.) The Discrete Fourier Transform is the same thing except it handles time and frequency samples rather than a continuous function.However, the Fourier Transfer (and therefore the DFT and FFT) work with complex numbers; which means twice as much storage and three times as many calculations than if real numbers were used. (Some of this can be reduced, but it adds complexity to the algorithm.) One alternative is the Hartley Transform and the Discrete Hartley (nee Bracewell) Transform which uses only real numbers.The DHT is very simple: Y[m] = SUM n=0..N-1 X[n] * cas(n*m*2pi/N) / sqrt(N)where N is the number of samples, cas(z) = cos(z) + sin(z), m = 0..N-1, and cas(a) is in radiansThe DHT is also it's own inverse (especially if expressed with the 1/sqrt(N) scaling factor), which is very nice as you can take the output and run it back through the same routine to get back to the original input. It's also trivial to get the Fourier sequence given the Hartley sequence:F[x] = ( H[x] + H[N-x] ) / 2 + ( H[x] - H[N-x] ) / 2iAnd just like the FFT, it's possible to calculate the DHT in O(N*log2(N)) operations. (Where N=2^P) This is from United States patent number 4,646,256, by Ronald N. Bracewell, assigned to Stanford University, and released to the public domain.

swap each element with the bitwise reversed address elementi.e. for N=16 element #5 (%0101) gets swapped with element #10 (%1010).for stage 1 to log2N  width = 2 ^ stage  span = width / 2  for i = 0 to span-1    cas = cos( i/width * 2pi ) + sin( i/width * 2pi )    for j = 0 to N-1 step width      left = sample[i+j]      right = sample[i+j+span] * cas      sample[i+j] = left + right      sample[i+j+span] = left - right    next j  next inext stage

And that's it! (Well, it's missing the scaling factor.) It should be noted stage = 1 and stage = 2 are trivial since cas == 1. And those who read the patent will discover that my algorithm looks nothing like the patent. This is partially because the patent is for an "implementation", and also because my algorithm minimizes the number of sin/cos functions.This works for N=4 but fails on N=8. I'm re-reading the patent and I think I see where I went wrong.