/* * * These functions can be assembled into various kinds of factored * discrete Fourier transforms. * * 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 "complex.h" #include "common.h" #include "bitrev.h" #include "fft.h" #define fftbitrev(out,in,q) bitrevd(out,in,q,sizeof(complex)) #define fftbrinpl(x,q) bitrevi(x,q,sizeof(complex)) /**************************************************************** * This function successively applies the sparse matrices * F_{q-1}, F_{q-2}, ..., F_1, and finally F_0 to the complex * input vector `f[]', transforming it in place. */ extern void fftproduct( /* Apply sparse matrix product. */ complex *f, /* Input and output vector. */ int q, /* Length of `f[]' is N=1<=0; k--) { N1 = N>>k; /* block size */ M = N1>>1; /* butterfly size */ /* Each F_k has 2^{k} blocks E_{2M}, where 2M=2^{q-k}, * and the matrix E_M is defined by: * I_M W_M * I_M -W_M */ for( b=0; bRe; /* Mult. f[M] by +/-W[0]= +/-1 */ tmp.Im = fptr2->Im; fptr2->Re = fptr1->Re - tmp.Re; fptr2->Im = fptr1->Im - tmp.Im; fptr1->Re += tmp.Re; fptr1->Im += tmp.Im; /* Do the cosine/sine terms */ Wptr = W; for( j=1; jRe = fptr1->Re - tmp.Re; fptr2->Im = fptr1->Im - tmp.Im; fptr1->Re += tmp.Re; fptr1->Im += tmp.Im; } } } return; } /**************************************************************** * Return a complex vector `Omega[]' of length `|M|' containing * the values: * W[n] = exp(PI*i*n/M), n = 0, 1,...,|M|-1. * If M<0, then the returned vector is the complex conjugate * of Omega(|M|). */ extern complex * fftomega( /* Return exp(-PI*i*n/M), */ int M) /* for n=0,1,2,...,|M|-1. */ { complex *W; double factor, theta; int n; factor = -PI/(double)M; M = absval(M); W = (complex *)malloc(M*sizeof(complex)); theta = 0.0; for(n=0; n= 0 ), or * - the inverse DFT of `f[]' * (if q < 0). * * External functions called: * calloc(),free() From ; allocate output vector. * fftbitrev() Permute `f[]' to `fhat[]' by bit-reversal. * fftnormal() Divide an `N'-component vector by sqrt(N). * fftomega() Return a vector of complex exponentials. * fftproduct() Apply the sparse matrices F_{q-1}...F_0. */ extern complex * dft( /* Allocate, assign and return */ const complex *f, /* a complex vector, the (1<