fftComplex

Computes the complex discrete Fourier transform of a complex sequence.

Synopsis

fftComplex (p)

Required Arguments

complex p[] (Input)
Array with n components containing the periodic sequence.

Return Value

If no optional arguments are used, fftComplex returns the transformed sequence. If no value can be computed, then None is returned.

Optional Arguments

backward
Compute the backward transform. If backward is used, the return value of the function is the backward transformed sequence.
params, float[] (Input)
Vector returned by a previous call to fftComplexInit. If fftComplex is used repeatedly with the same value of n, then it is more efficient to compute these parameters only once.

Description

The function fftComplex 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, fftComplex computes the forward transform below.

\[q_j = \sum_{m=0}^{n-1} p_m e^{-2 \pi imj / n}\]

Note that we can invert the Fourier transform as follows below.

\[p_m = \tfrac{1}{n} \sum_{j=0}^{n-1} q_j e^{2 \pi ijm/n}\]

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 backward is selected, then the following computation is performed.

\[q_j = \sum_{m=0}^{n-1} p_m e^{2 \pi imj/n}\]

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 \(p_2 = np\).

q  = fftComplex(p)
p2 = fftComplex(q, backward)

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.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.fftComplex import fftComplex

n = 7
two_pi = 2 * constant("pi")
p = empty(n, dtype=complex)

# Fill p with a pure exponential signal
for k in range(0, n):
    p[k] = exp(complex(0.0, k * two_pi / n))

q = fftComplex(p)

print("        index    p.re      p.im      q.re      q.im")
for k in range(0, n):
    print("%11d%10.2f%10.2f%10.2f%10.2f"
          % (k, p[k].real, p[k].imag, q[k].real, q[k].imag))

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.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.fftComplex import fftComplex

n = 7
two_pi = 2 * constant("pi")
p = empty(n, dtype=complex)

# Fill p with an increasing signal
for k in range(0, n):
    p[k] = complex(k, 0)

q = fftComplex(p)
pp = fftComplex(q, backward=True)

print("        index    p.re      p.im     pp.re     pp.im")
for k in range(0, n):
    print("%11d%10.2f%10.2f%10.2f%10.2f"
          % (k, p[k].real, p[k].imag, pp[k].real, pp[k].imag))

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