fftReal¶

Computes the real discrete Fourier transform of a real sequence.

Synopsis¶

fftReal (p)

Required Arguments¶

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

Return Value¶

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 fftRealInit. If fftReal is used repeatedly with the same value of n, then it is more efficient to compute these parameters only once.

Description¶

The function fftReal 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 real FFT in FFTPACK, which was developed by Paul Swarztrauber at the National Center for Atmospheric Research.

By default, fftReal computes the forward transform. If n is even, then the forward transform is

\[\begin{split}\begin{array}{l} q_{2m-1} = \displaystyle\sum_{k=0}^{n-1} p_k\cos\tfrac{2{\pi}km}{n}\phantom{...}m = 1, \ldots, n / 2 \\ q_{2m} = - \displaystyle\sum_{k=0}^{n-1} p_k\sin\tfrac{2{\pi}km}{n}\phantom{...}m = 1, \ldots, n/2 - 1 \\ q_0 = \displaystyle\sum_{k=0}^{n-1} p_k \\ \end{array}\end{split}\]

If n is odd, \(q_m\) is defined as above for m from 1 to \((n − 1)/2\).

Let f be a real valued function of time. Suppose we sample f at n equally spaced time intervals of length Δ seconds starting at time \(t_0\). That is, we have

\[p_i:= f(t_0 + iΔ) i = 0, 1, …, n − 1\]

We will assume that n is odd for the remainder of this discussion. The function fftReal treats this sequence as if it were periodic of period n. In particular, it assumes that \(f(t_0) = f(t_0 + n * \Delta)\). Hence, the period of the function is assumed to be \(T = n * \Delta\). We can invert the above transform for p as follows:

\[p_m = \frac{1}{n} \left[ q_0 + 2\sum_{k=0}^{(n-3)/2} q_{2k+1}\cos\frac{2\pi(k+1)m}{n} - 2\sum_{k=0}^{(n-3)/2} q_{2k+2}\sin\frac{2\pi(k+1)m}{n} \right]\]

This formula is very revealing. It can be interpreted in the following manner. The coefficients q produced by fftReal determine an interpolating trigonometric polynomial to the data. That is, if we define

\[\begin{split}\begin{aligned} g(t) &= \frac{1}{n} \left[ q_0 + 2 \sum_{k=0}^{(n-3)/2} q_{2k+1}\cos\frac{2\pi(k+1)(t-t_0)}{n\mathit{\Delta}}-2\sum_{k=0}^{(n-3)/2}q_{2k+2}\sin\frac{2\pi(k+1)(t-t_0)}{n\mathit{\Delta}} \right] \\ &= \frac{1}{n} \left[ q_0 + 2 \sum_{k=0}^{(n-3)/2} q_{2k+1}\cos\frac{2\pi(k+1)(t-t_0)}{T}-2\sum_{k=0}^{(n-3)/2}q_{2k+2}\sin\frac{2\pi(k+1)(t-t_0)}{T} \right] \\ \end{aligned}\end{split}\]

then we have

\[f \left(t_0 + i\mathit{\Delta}\right) = g\left(t_0 + i\mathit{\Delta}\right)\]

Now suppose we want to discover the dominant frequencies, forming the vector P of length (n + 1)/2 as follows:

\[\begin{split}\begin{array}{l} P_0 := \left|q_0\right| \\ P_k := \sqrt{q_{2k-1}^2 + q_{2k}^2} \phantom{...} k = 1,2, \ldots, (n-1)/2 \end{array}\end{split}\]

These numbers correspond to the energy in the spectrum of the signal. In particular, \(P_k\) corresponds to the energy level at frequency

\[\frac{k}{T} = \frac{k}{n\mathit{\Delta}} \phantom{...} k=0,1,\ldots,\frac{n-1}{2}\]

Furthermore, note that there are only (n + 1)/2 ≈ T/(2Δ) resolvable frequencies when n observations are taken. This is related to the Nyquist phenomenon, which is induced by discrete sampling of a continuous signal. Similar relations hold for the case when n is even.

If the optional argument backward is specified, then the backward transform is computed. If n is even, then the backward transform is

\[q_m = p_0 + (-1)^mp_{n-1} + 2\sum_{k=0}^{n/2-2}p_{2k+1}\cos\frac{2\pi(k+1)m}{n} - 2\sum_{k=0}^{n/2-2}p_{2k+2}\sin\frac{2\pi(k+1)m}{n}\]

If n is odd,

\[q_m = p_0 + 2 \sum_{k=0}^{(n-3)/2} p_{2k+1} \cos \frac{2\pi(k+1)m}{n} - 2 \sum_{k=0}^{(n-3)/2} p_{2k+2} \sin \frac{2\pi(k+1)m}{n}\]

The backward Fourier transform is the unnormalized inverse of the forward Fourier transform.

Examples¶

Example 1¶

In this example, a pure cosine wave is used as a data vector, and its Fourier series is recovered. The Fourier series 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.fftReal import fftReal

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

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

q = fftReal(p)

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

Output¶

        index      p        q
          0      1.00     -0.00
          1      0.62      3.50
          2     -0.22     -0.00
          3     -0.90      0.00
          4     -0.90     -0.00
          5     -0.22     -0.00
          6      0.62     -0.00

Example 2¶

This example computes the Fourier transform of the vector x, where \(x_j = (-1)^j\) for \(j = 0\) to \(n − 1\). The backward transform of this vector is now computed by using the optional argument backward. Note that \(s = nx\), that is, \(s_j = (-1)^j n\), for \(j = 0\) to \(n − 1\).

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

n = 7
two_pi = 2 * constant("pi")
x = empty(n)

# Fill data vector
x[0] = 1.0
for k in range(1, n):
    x[k] = -x[k - 1]

# Compute the forward transform of x
q = fftReal(x)

# Compute the backward transform of x
s = fftReal(q, backward=True)
print("        index      x        q        s")
for k in range(0, n):
    print("%11d%10.2f%10.2f%10.2f" % (k, x[k], q[k], s[k]))

Output¶

        index      x        q        s
          0      1.00      1.00      7.00
          1     -1.00      1.00     -7.00
          2      1.00      0.48      7.00
          3     -1.00      1.00     -7.00
          4      1.00      1.25      7.00
          5     -1.00      1.00     -7.00
          6      1.00      4.38      7.00