FFT and Fast Fourier Sine Transform

CCATSL provides two routines for calculating transforms of discrete sequences, both based on the Fast Fourier Transform. FftCL implements the basic Fast Fourier algorithm, while FftSinCL computes a Sine transform using FftCL.

Fast Fourier Transform, FftCL

FftCL computes the Discrete Fourier Transform of a sequence $X_0,X_1,\ldots{,}\,X_{N-1}$ of complex numbers, using the Fast Fourier algorithm. The algorithm is an efficient way of calculating the transformed sequence:

$$\m \widetilde X_s = \sum_{r=0}^{N-1} e^{-2\pi i s r /N}X_r$$

when $N=2^m$ for some $m$. The original sequence can be recovered from the transformed sequence $\widetilde X_s$ by

$$\m X_r = \frac{1}{N}\sum_{s=0}^{N-1} e^{2\pi i r s /N}\widetilde X_s$$

which is essentially just another application of the FFT algorithm (followed by division by $N$), but with a change of sign in the exponential.

	void FftCL (	double *x,
		int m,
		double sign);

	x		Array holding the sequence $X_0,\ldots{,}\,X_{N-1}$. The array x points to must be an double [0..N-1][0..1] where x[i][0] is the real part of $X_i$ and x[i][1] is the imaginary part. FftCL replaces these values with the transformed sequence, $\widetilde X_0,\ldots{,}\,\widetilde X_{N-1}$. See below for the an example of how to setup a suitable array.
	m		Specifies the number of terms in the sequence; there must be $N$ = $2^n$ of them.
	sign		The sign of the exponential in the series above (-1 for the normal transform, 1 for the inverse). Note that when computing the inverse fast Fourier transform by using FftCL with sign=1 you will need to manually divide each resulting $X_i$ by $N$ (see the example below).

Here is an example using FftCL:

/* /examples/chapter2/fftex.c */ #include <catam.h> int MainCL(void) { int r=3; int s; double x[8][2]; for (s=0; s<8; s++) { /* Setup X to be the rth Fourier mode. Here the rth Fourier mode, F_s is given by F_s = (1/8)exp(2 PI i R s / 8) (where i=sqrt(-1)) which we split into its real and imaginary parts. */ x[s][0]=cos(2*M_PI*r*s/8.0)/8.0; x[s][1]=sin(2*M_PI*r*s/8.0)/8.0; printf("x[%i]=(%f + i %f)",s,x[s][0],x[s][1]); } printf(""); FftCL(x,3,-1.0); /* transform */ for (s=0; s<8; s++) printf("x[%i]=(%f +i %f)",s,x[s][0],x[s][1]); FftCL(x,3,1.0); /* transform back*/ printf(""); for (s=0; s<8; s++) printf("x[%i]=(%f +i %f)",s,x[s][0]/8.0,x[s][1]/8.0); return 0; }

Fast Fourier Sine Transform, FftSinCL

FftSinCL computes the sine transform of a sequence $X_0,X_1,\ldots{,}\,X_{N-1},X_N$ of real numbers with $X_0$ = $X_N$ = 0, where $N$ = $2^m$ for some $m$:

$$\m \widetilde X_s=\sum_{r=0}^{N-1} \sin(\pi rs/N)X_r.$$

It can also compute the inverse transform. (Note that unlike FftCL, when FftSinCL performs the inverse transform, no subsequent division by $N$ is necessary.)

	void FftSinCL (	int m,
		double *x,
		double dirn);

	m		Specifies the number of terms in the sequence; since $N=2^m$, there will be $2^m+1$ terms.
	x		Array holding the sequence $X_0,\ldots{,}\,X_N$. FftSinCL will replace these values with the transformed sequence.
	dirn		Indicates whether the normal sine transform (dirn=1) or the inverse transform (dirn=-1) is required.

A typical use of FftSinCL might look like:

/* /examples/chapter2/sinfftex.c */ .. { double x[256]; /* setup the array x */ .. FftSinCL(8,x,1.0); }

A demonstration of FftSinCL can be found in b08poisn.c.