The FLIBSCI.LIB Fourier and Wavelet Transform Software

Object modules to perform wavelet and fourier transforms are to be found in the library flibsci.lib. This library is compatible with the Microsoft Visual C/C++ compiler and linker.

The Fourier Transform Software

Object modules have been prepared to execute a discrete complex fourier transform. Also supplied are routines to calculate inverse fourier transforms, convolution, and correlation.

The fast Fourier transform routine has the entry point named FFT. There are no subroutine arguments and no function value is returned.

The preparatory subroutine FFTPR is called first, and only once, that fills arrays with the appropriate sine and cosine values needed by the FFT routine. These values do not change from calculation to calculation, so they need be calculated only once.

The inverse fourier transform, the auto correlation, cross correlation and convolution of 8, 16, 32, 64, 128, 256, 512, and 1024, 2048 and 4096 point complex data sets all call the routine FFT. The routine FFTPR and the other auxilliary routines are to be found in the source modules FFTETC.FOR (and FFTETC.C, if desired).

The inverse fourier transform is calculated by performing the fourier transform of the complex conjugate of the function Z(f) and then taking the complex conjugate of the result (by changing the sign of the imaginary part) and then dividing by the number of sample points, N. This process is carried out by the routine IFFT which is supplied in source module form in Fortran and C in the files FFTETC.FOR and .C, respectively.

The routines given in FFTETC.FOR to calculate convolution (CONV), autocorrelation (AUTOCORR) and cross correlation (CROSS) are all written in Fortran (or C if FFTETC.C is used). Each of these calls IFFT and FFT to do the work. The routine IFFT calculates the inverse fourier transform and it in turn calls FFT, which is the only routine written in assembly language.

Remember that the discrete transform yields correct results for continuous functions only as long as two important restrictions are obeyed. These are:

1. The function to be transformed must be periodic with a period equal to the sampling interval.

2. The sampling rate must be more than twice the maximum frequency content of the function to be transformed.

Otherwise, the user is destined to learn more about aliasing, the Gibbs phenomenon and other things than he probably ever wanted to know.

The wavelet transformation has no such restrictions!

The routine SUBROUTINE FFTPR(NU), included with the FFT software in the library FLIBSCI.LIB, must be called before any of the Fourier transform routines, that are described below, may be used. Furthermore, it must be called again whenever the value of the variable NU is changed. Unless the value of NU is changed in a program, this routine need only be called once; however, it does no harm to call this routine any number of times.

The number of points utilized by each of the routines is given by specifying the value of NU which is the LOG to the base 2 of the number of sampled data points, NFFT. That is, NFFT=2**NU.

Each of the routines uses certain data structures, irrespective of the number of data points included in the transform calculation. Clearly, the number of data points is a power of 2. The maximum number of points is 4096; that is, NU is less than or equal to 12.

The data structures (COMMON blocks) that are used in each of the routines is as follows:

 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMPLEX*16 XC
      COMMON/FFTT/XC(4096)
      COMMON/NFFT/NFFT

The values of the sampled data points are to be inserted into and read from the array XC.

The routine SUBROUTINE FFT calculates the (forward) fourier transform on the array of values given in the complex array XC. The results are returned in the array XC, therefore the initial values in the array XC are destroyed in this process.

The routine SUBROUTINE IFFT calculates the (inverse) fourier transform on the array of values given in the complex array XC. The results are returned in the array XC, therefore the initial values in the array XC are destroyed in this process.

The routine SUBROUTINE CONV(G) calculates the convolution of two complex functions, XC and G. The result is placed in XC. G is an array of COMPLEX*16 values.

The routine SUBROUTINE AUTOCOR calculates the auto-correlation function of the complex function in XC. The result is placed in XC.

The routine SUBROUTINE CROSS(G) calculates the cross-correlation function of the complex functions in XC and G. The result is placed XC. G is an array of COMPLEX*16 values.

The routine SUBROUTINE HARTLEY calculates the Hartley transform of the real function given in the array XC. The result is placed in the array XC. The array XC is destroyed in this process.

The routine SUBROUTINE IHARTLEY calculates the inverse Hartley transform of the real function given in the array XC. The result is placed in the array XC. The origional values of the array XC is destroyed in this process.

The Forward and Inverse Wavelet Transformation Software

Software routines have been implemented to perform the forward and inverse transformations for the Haar (M = 2) and the D4 wavelets (M = 4). These have been written to handle 2, 4, 8, 16, 32, 64, 128, 256, 512 and 1024 data points.

The forward and inverse wavelet transform routines use the data structures given here:

 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/WAVELET/F(4096),A(8192),B(4096),C(0:6)

These data structures (COMMON blocks) are used irrespective of the number of data points ,N, that are employed in performing wavelet transformations, the type of wavelet used, and the number of data points employed.

The value of N, the number of points used, is calculated by the wavelet transform routines from the value of JX that is given as one of the subroutine arguments. The value of JX is the LOG to the base 2 of the number of points. That is, N=2**JX. This forces the number of points to be a power of 2.

The value of JX must be less than of equal to 12.

The value of MX must be either 2, for Haar wavelets, or 4 for D4 wavelets.

A COMMON block. /WAVETC/, is defined within the wavelet transform routines. It is for the private use of the FLIBSCI routines. The values may be read but not modified by anything outside the wavelet transform routines. Its structure is given as follows, for reference only:

 
      COMMON/WAVEETC/N,N2,NT,JT,M

The routine SUBROUTINE WAVE(MX,JX) carries out the (forward) wavelet transformation on a real function given in the array F held in the COMMON block /WAVELET/. The results are placed in the arrays A and B.

The way to cause a forward wavelet transformation to be performed is to use

 
      CALL WAVE(M,JT)

where N = 2**JT; i.e., for JT = 3 then N = 8, for JT = 4 then N = 16 , etc. The array, F, must have been filled with the data to be transformed before the call is made. The transform is performed and the A and B arrays are filled with the "answers".

The routine SUBROUTINE IWAVE(MX,JX) carries out the inverse wavelet transformation on a real function given in the array B held in the COMMON block /WAVELET/. The results are placed in the arrays A and F.

The inverse transformation is performed by using

 
     CALL IWAVE(M,JT)

The B array must have been filled with data to be back transformed. The arrays F and A are calculated by IWAVE. Again, M is the number of wavelet coefficients and the number of sampled data points is equal to 2**JT.