Hartley Transform

L = log2(N)for i = 0 to N-1  j = bitreverse( i, L )  if j > i then    t = x[i]    x[i] = x[j]    x[j] = t  endifnext ifor i = 0 to N-1 step 4  t0 = x[i]  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[i]    t1 = x[i+q]    t2 = x[i+h]    t2 = x[i+h+q]    x[i] = 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