/* * * Functions to perform the discrete Hartley transform. * * Copyright (C) 1991--94 Wickerhauser Consulting. All Rights Reserved. * May be freely copied for personal, noncommercial use by owners of * ``Adapted Wavelet Analysis from Theory to Software'' ISBN 1-56881-041-5 * by Mladen Victor Wickerhauser [AK Peters, Ltd., Wellesley, Mass., 1994] */ #include #include #include #include "real.h" #include "bitrev.h" #include "hartley.h" #define SQH (0.707106781186547) /* Square root of 1/2 */ #define PI (3.1415926535897932) #define dhtbitrev(out,in,q) bitrevd(out,in,q,sizeof(real)); #define dhtbrinpl(x,q) bitrevi(x,q,sizeof(real)); /**************************************************************** * This function applies the sparse matrices G_{q-1}, G_{q-2}, * ... G_1, and finally G_0 to the real input vector `f[]' in * place. We assume that q>2. */ extern void dhtproduct( /* Apply sparse matrix product. */ real *f, /* Input and output vector. */ int q, /* Length of `f[]' is N=1<2); /* Test our assumption. */ N= 1<=0; k--) { N1 = N>>k; /* size of blocks in G_k */ M = N1>>1; /* butterfly size in G_k */ M2 = M>>1; /* midpoint of butterfly */ /* Each G_k is the block diagonal direct sum of * 2^{k} blocks A_{2M}, each of size N1 x N1, where * N1 = 2M = 2^{q-k} and the matrix A_{2M} is * defined by: * I_M B_M * I_M -B_M * * The matrix B_M is the "butterfly" below: * * 1 0 . . . 0 * 0 c1 s1 * c2 s2 * . . * . cK sK * . 1 * . sK -cK * . . * s2 -c2 * 0 s1 -c1 * * for K=(M/2)-1, ck=cos(PI*k/M), sk=sin(PI*k/M). * */ for( b=0; b>1); n++) { theta += factor; c[n] = (real)cos(theta); s[n] = (real)sin(theta); } return; } /**************************************************************** * Normalize `f[]', a vector of length `N', by dividing each * component by sqrt(N). */ extern void dhtnormal( real *f, /* Multiply `f[n]' for n=0,1, */ int N ) /* ...,N by `1.0/sqrt(N).' */ { double norm; int n; norm = sqrt(1.0/(double)N); for(n=0; n; to allocate arrays. * dhtbitrev() Permute `f[]' to `fhat[]' by bit-reversal. * dhtnormal() Divide an `N'-component vector by sqrt(N). * dhtcossin() Fill arrays with sines and cosines.. * dhtproduct() Apply the sparse matrices F_{q-1}...F_0. */ extern real * dht( /* Allocate, assign and return */ const real *f, /* a real vector, the (1<=0); /* Accept only positive `q'. */ N = 1<