AtariAge Forums: AtariAge Forums -> Hartley Transform

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 radians

The 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] ) / 2i

And 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 element
i.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 i
next 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.