fft_complex
Computes the complex discrete Fourier transform of a complex sequence.
Synopsis
#include <imsl.h>
f_complex *imsl_c_fft_complex (int n, f_complex p[], …, 0)
The type d_complex function is imsl_z_fft_complex.
Required Arguments
int n (Input)
Length of the sequence to be transformed.
f_complex p[] (Input)
Array with n components containing the periodic sequence.
Return Value
If no optional arguments are used, imsl_c_fft_complex returns a pointer to the transformed sequence. To release this space, use imsl_free. If no value can be computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
f_complex *imsl_c_fft_complex (int n, f_complex p[],
IMSL_BACKWARD,
IMSL_PARAMS, float params[],
IMSL_RETURN_USER, f_complex q[],
0)
Optional Arguments
IMSL_BACKWARD
Compute the backward transform. If IMSL_BACKWARD is used, the return value of the function is the backward transformed sequence.
IMSL_PARAMS, float params[] (Input)
Pointer returned by a previous call to imsl_c_fft_complex_init. If imsl_c_fft_complex is used repeatedly with the same value of n, then it is more efficient to compute these parameters only once.
IMSL_RETURN_USER, f_complex q[] (Output)
Store the result in the user-provided space pointed to by q. Therefore, no storage is allocated for the solution, and imsl_c_fft_complex returns q. The array q must be of length at least n.
Description
The function imsl_c_fft_complex computes the discrete Fourier transform of a real vector of size n. It uses the system’s high performance library for the computation, if available. Otherwise, the method used is a variant of the Cooley-Tukey algorithm, which is most efficient when n is a product of small prime factors. If n satisfies this condition, then the computational effort is proportional to n log n. The Cooley-Tukey algorithm is based on the complex FFT in FFTPACK, which was developed by Paul Swarztrauber at the National Center for Atmospheric Research.
By default, imsl_c_fft_complex computes the forward transform below.
Note that we can invert the Fourier transform as follows below.
This formula reveals the fact that, after properly normalizing the Fourier coefficients, you have the coefficients for a trigonometric interpolating polynomial to the data.
If the option IMSL_BACKWARD is selected, then the following computation is performed.
Furthermore, the relation between the forward and backward transforms is that they are unnormalized inverses of each other. That is, the following code fragment begins with a vector p and concludes with a vector p2 = np.
q = imsl_c_fft_complex(n, p, 0);
p2 = imsl_c_fft_complex(n, q, IMSL_BACKWARD, 0);
Examples
Example 1
This example inputs a pure exponential data vector and recovers its Fourier series, which is a vector with all components zero except at the appropriate frequency where it has an n.
#include <imsl.h>
#include <stdio.h>
int main()
{
int k, n = 7;
float two_pi = 2*imsl_f_constant("pi", 0);
f_complex p[7], *q, z;
/* Fill p with a pure exponential signal */
for (k = 0; k < n; k++) {
z.re = 0.;
z.im = k*two_pi/n;
p[k] = imsl_c_exp(z);
}
q = imsl_c_fft_complex (n, p, 0);
printf(" index p.re p.im q.re q.im\n");
for (k = 0; k < n; k++)
printf("%11d%10.2f%10.2f%10.2f%10.2f\n", k, p[k].re, p[k].im,
q[k].re, q[k].im);
}
Output
index p.re p.im q.re q.im
0 1.00 0.00 0.00 0.00
1 0.62 0.78 7.00 0.00
2 -0.22 0.97 -0.00 0.00
3 -0.90 0.43 -0.00 0.00
4 -0.90 -0.43 0.00 -0.00
5 -0.22 -0.97 0.00 -0.00
6 0.62 -0.78 0.00 0.00
Example 2
The backward transform is used to recover the original sequence. Notice that the forward transform followed by the backward transform multiplies the entries in the original sequence by the length of the sequence.
#include <imsl.h>
#include <stdio.h>
int main()
{
int k, n = 7;
float two_pi = 2*imsl_f_constant("pi", 0);
f_complex p[7], *q, *pp;
/* Fill p with an increasing signal */
for (k = 0; k < n; k++) {
p[k].re = (float) k;
p[k].im = 0.;
}
q = imsl_c_fft_complex (n, p, 0);
pp = imsl_c_fft_complex (n, q,
IMSL_BACKWARD,
0);
printf(" index p.re p.im pp.re pp.im \n");
for (k = 0; k < n; k++)
printf("%11d%10.2f%10.2f%10.2f%10.2f\n", k, p[k].re, p[k].im,
pp[k].re , pp[k].im);
}
Output
index p.re p.im pp.re pp.im
0 0.00 0.00 0.00 0.00
1 1.00 0.00 7.00 0.00
2 2.00 0.00 14.00 0.00
3 3.00 0.00 21.00 0.00
4 4.00 0.00 28.00 0.00
5 5.00 0.00 35.00 0.00
6 6.00 0.00 42.00 0.00
