fft_real
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 imsl_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. 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_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. 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, 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
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.
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>
int main()
{
int k, n = 7;
float two_pi = 2*imsl_f_constant("pi", 0);
float p[7], *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>
int main()
{
int k, n = 7;
float *q, *s, x[7];
/* 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
