EricBall's Tech Projects - Hartley Transform

[RSS Bot]

Sample code (unoptimized). Comment if you want more details.

 

L = log2(N)for i = 0 to N-1 j = bitreverse( i, L ) if j > i then t = x x = x[j] x[j] = t endifnext ifor i = 0 to N-1 step 4 t0 = x t2 = x[i+1] t1 = x[i+2] t3 = x[i+3] x[i+0]=(t0+t1+t2+t3)/2 x[i+1]=(t0+t1-t2-t3)/2 x[i+2]=(t0-t1+t2-t3)/2 x[i+3]=(t0-t1-t2+t3)/2next ifor o = 3 to L s = 1 << o h = s / 2 q = s / 4 for i = 0 to N-1 step s t0 = x t1 = x[i+q] t2 = x[i+h] t2 = x[i+h+q] x = t0 + t2 x[i+q] = t1 + t3 x[i+h] = t0 - t2 x[i+h+q] = t1 - t3 for j = 1 to h-1 t0 = x[i+j] / sqrt(2) t1 = x[i+h-j] / sqrt(2) t2 = x[i+h+j] / sqrt(2) t2 = x[i+s-j] / sqrt(2) x[i+j] = t0 + t1 * cos(2*pi()*j/N) + t3 * sin(2*pi()*j/N) x[i+h-j] = t2 + t1 * sin(2*pi()*j/N) - t3 * cos(2*pi()*j/N) x[i+h+j] = t0 - t1 * cos(2*pi()*j/N) - t3 * sin(2*pi()*j/N) x[i+s-j] = t2 - t1 * sin(2*pi()*j/N) + t3 * cos(2*pi()*j/N) next j next inext o

 

http://www.atariage.com/forums/index.php?app=blog&blogid=7&showentry=6654