fft_2d_complex
Computes the complex discrete two-dimensional Fourier transform of a complex two-dimensional array.
Synopsis
#include <imsl.h>
f_complex *imsl_c_fft_2d_complex (int n, int m, f_complex p[], …, 0)
The type d_complex function is imsl_z_fft_2d_complex.
Required Arguments
int n (Input)
Number of rows in the two-dimensional transform.
Number of rows in the two-dimensional transform.
int m (Input)
Number of columns in the two-dimensional transform.
Number of columns in the two-dimensional transform.
f_complex p[] (Input)
Two-dimensional array of size n × m containing the sequence that is to be transformed.
Two-dimensional array of size n × m containing the sequence that is to be transformed.
Return Value
A pointer to the transformed array. 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_2d_complex (int n, int m, f_complex p[],
IMSL_P_COL_DIM, int p_col_dim,
IMSL_BACKWARD,
IMSL_RETURN_USER, f_complex q[],
IMSL_Q_COL_DIM, int q_col_dim,
0)
Optional Arguments
IMSL_P_COL_DIM, int p_col_dim (Input)
The column dimension of p.
Default: p_col_dim = m
The column dimension of p.
Default: p_col_dim = m
IMSL_BACKWARD
Compute the backward transform. If IMSL_BACKWARD is used, the return value of the function is the backward transformed sequence.
Compute the backward transform. If IMSL_BACKWARD is used, the return value of the function is the backward transformed sequence.
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_2d_complex returns q. The array must be of length at least n × m.
Store the result in the user-provided space pointed to by q. Therefore, no storage is allocated for the solution, and imsl_c_fft_2d_complex returns q. The array must be of length at least n × m.
IMSL_Q_COL_DIM, int q_col_dim (Input)
The column dimension of q.
Default: q_col_dim = m
The column dimension of q.
Default: q_col_dim = m
Description
The function imsl_c_fft_2d_complex computes the discrete Fourier transform of a two-dimensional complex array of size n × m. 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 both n and m are a product of small prime factors. If n and m satisfy this condition, then the computational effort is proportional to nm log nm. 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_2d_complex computes the forward transform below.
Note that we can invert the Fourier transform as follows.
This formula reveals the fact that, after properly normalizing the Fourier coefficients, you have the coefficients for a trigonometric interpolating polynomial to the data. The function imsl_c_fft_2d_complex is based on the complex FFT in FFTPACK, which was developed by Paul Swarztrauber at the National Center for Atmospheric Research.
If the option IMSL_BACKWARD is selected, then the following computation is performed.
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 = nmp.
q = imsl_c_fft_2d_complex(n, m, p, 0);
p2 = imsl_c_fft_2d_complex(n, m, q, IMSL_BACKWARD, 0);
Examples
Example 1
This example computes the Fourier transform of the pure frequency input for a 5 × 4 array
for 0 ≤ n ≤ 4 and 0 ≤ m ≤ 3. The result, , has all zeros except in the [2][3] position.
#include <imsl.h>
int main()
{
int s, t, n = 5, m =4;
float two_pi = 2*imsl_f_constant("pi", 0);
f_complex p[5][4], *q, z, w;
/* Fill p with a pure exponential signal */
for (s = 0; s < n; s++) {
z.re = 0.;
z.im = s*two_pi*2./n;
for(t =0; t < m; t++){
w.re = 0.;
w.im = t*two_pi*3./m;
p[s][t] = imsl_c_mul(imsl_c_exp(z),imsl_c_exp(w));
}
}
q = imsl_c_fft_2d_complex (n, m, (f_complex*)p,
0);
/* Write the input */
imsl_c_write_matrix ("The input matrix is ", 5, 4, (f_complex*)p,
IMSL_ROW_NUMBER_ZERO,
IMSL_COL_NUMBER_ZERO,
0);
imsl_c_write_matrix ("The output matrix is ", 5, 4, q,
IMSL_ROW_NUMBER_ZERO,
IMSL_COL_NUMBER_ZERO,
0);
}
Output
The input matrix is
0 1
0 ( 1.000, 0.000) ( 0.000, -1.000)
1 ( -0.809, 0.588) ( 0.588, 0.809)
2 ( 0.309, -0.951) ( -0.951, -0.309)
3 ( 0.309, 0.951) ( 0.951, -0.309)
4 ( -0.809, -0.588) ( -0.588, 0.809)
3 ( -0.309, -0.951) ( -0.951, 0.309)
4 ( 0.809, 0.588) ( 0.588, -0.809)
2 3
0 ( -1.000, -0.000) ( -0.000, 1.000)
1 ( 0.809, -0.588) ( -0.588, -0.809)
2 ( -0.309, 0.951) ( 0.951, 0.309)
The output matrix is
0 1
0 ( -0, -0) ( 0, -0)
1 ( 0, 0) ( 0, -0)
2 ( -0, -0) ( 0, -0)
3 ( 0, 0) ( 0, -0)
4 ( -0, -0) ( 0, -0)
2 3
0 ( 0, -0) ( 0, -0)
1 ( -0, 0) ( 0, -0)
2 ( 0, -0) ( 20, 0)
3 ( -0, 0) ( -0, -0)
4 ( 0, -0) ( -0, -0)
Example 2
This example uses the backward transform to recover the original sequence. Notice that the forward transform followed by the backward transform multiplies the entries in the original sequence by the product of the lengths of the two dimensions.
#include <imsl.h>
#include <stdio.h>
int main()
{
int s, t, n = 5, m =4;
f_complex p[5][4], *q, *p2;
/* Fill p with a pure exponential signal */
for (s = 0; s < n; s++) {
for(t =0; t < m; t++){
p[s][t].re = s + 5*t;
p[s][t].im = 0.;
}
} /* Forward transform */
q = imsl_c_fft_2d_complex (n, m, (f_complex*)p, 0);
/* Backward transform */
p2 = imsl_c_fft_2d_complex (n, m, q,
IMSL_BACKWARD, 0);
/* Write the input */
imsl_c_write_matrix ("The input matrix is ", 5, 4, (f_complex*)p,
IMSL_ROW_NUMBER_ZERO,
IMSL_COL_NUMBER_ZERO, 0);
imsl_c_write_matrix ("The output matrix is ", 5, 4, p2,
IMSL_ROW_NUMBER_ZERO,
IMSL_COL_NUMBER_ZERO, 0);
}
Output
The input matrix is
0 1
0 ( 0, 0) ( 5, 0)
1 ( 1, 0) ( 6, 0)
2 ( 2, 0) ( 7, 0)
3 ( 3, 0) ( 8, 0)
4 ( 4, 0) ( 9, 0)
2 3
0 ( 10, 0) ( 15, 0)
1 ( 11, 0) ( 16, 0)
2 ( 12, 0) ( 17, 0)
3 ( 13, 0) ( 18, 0)
4 ( 14, 0) ( 19, 0)
The output matrix is
0 1
0 ( 0, 0) ( 100, 0)
1 ( 20, 0) ( 120, 0)
2 ( 40, 0) ( 140, 0)
3 ( 60, 0) ( 160, 0)
4 ( 80, 0) ( 180, 0)
2 3
0 ( 200, 0) ( 300, 0)
1 ( 220, 0) ( 320, 0)
2 ( 240, 0) ( 340, 0)
3 ( 260, 0) ( 360, 0)
4 ( 280, 0) ( 380, 0)
