Computes the real discrete Fourier transform of a real sequence.

Synopsis

#include <imsl.h>

float *imsl_f_fft_real (int n, float p[], ¼, 0)

The type double function is imsl_d_fft_real.

Required Arguments

int n   (Input)
Length of the sequence to be transformed.

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

Return Value

A pointer to the transformed sequence. To release this space, use free. If no value can be computed, then NULL is returned.

Synopsis with Optional Arguments

#include <imsl.h>

float *imsl_f_fft_real (int n, float p[],
IMSL_BACKWARD,
IMSL_PARAMS, float params[],
IMSL_RETURN_USER, float q[],
0)

Optional Arguments

IMSL_BACKWARD
Compute the backward transform and return a pointer to the (backward) transformed sequence.

IMSL_PARAMS, float params[]   (Input)
Pointer returned by a previous call to imsl_f_fft_real_init. If imsl_f_fft_real is used repeatedly with the same value of n, then it is more efficient to compute these parameters only once.

IMSL_RETURN_USER, float q[]   (Output)
Store the result in the user-provided space pointed to by q. Therefore, no storage is allocated for the solution, and imsl_f_fft_real returns q. The array q must be at least n long.

Description

The function imsl_f_fft_real computes the discrete Fourier transform of a real vector of size n. 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.

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

If n is odd, qm 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 t0. That is, we have

pi:= f(t0 + iΔ)    i = 0, 1, …, n − 1

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

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

then we have

f(t0 + iΔ) = g(t0 +iΔ)

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

These numbers correspond to the energy in the spectrum of the signal. In particular,
Pk corresponds to the energy level at frequency

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 IMSL_BACKWARD is specified, then the backward transform is computed. If n is even, then the backward transform is

If n is odd,

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

The function imsl_f_fft_real is based on the real FFT in FFTPACK, which was developed by Paul Swarztrauber at the National Center for Atmospheric Research.

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.

#include <imsl.h>
#include <math.h>
#include <stdio.h>

main()
{
    int         k, n = 7;
    float       two_pi = 2*imsl_f_constant("pi", 0);
    float       p[8], *q;
                                /* Fill q with a pure exponential signal */
    for (k = 0;  k < n;  k++)
        p[k] = cos(k*two_pi/n);

    q = imsl_f_fft_real (n, p, 0);

    printf("        index      p        q\n");
    for (k = 0;  k < n;  k++)
        printf("%11d%10.2f%10.2f\n", 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 xj = (−1)j for j = 0 to n − 1. The backward transform of this vector is now computed by using the optional argument IMSL_BACKWARD. Note that s = nx, that is,  sj = (−1)jn, for j = 0 to n − 1.

#include <imsl.h>
#include <stdio.h>

main()
{
    int         k, n = 7;
    float       *q, *s, x[8];
                                /* Fill data vector */
    x[0] = 1.0;
    for (k = 1;  k<n;  k++)
        x[k] = -x[k-1];
                                /* Compute the forward transform of x */
    q = imsl_f_fft_real (n, x, 0);
                               /* Compute the backward transform of x */
    s = imsl_f_fft_real (n, q,
                         IMSL_BACKWARD,
                         0);
    printf("        index      x        q        s\n");
    for (k = 0;  k < n;  k++)
        printf("%11d%10.2f%10.2f%10.2f\n", 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

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260