Introduction

The vDSP API provides Fourier transforms for transforming one-dimensional and two-dimensional data between the time domain and the frequency domain.

FFT Weights Arrays

To boost performance, vDSP functions that process frequency-domain data expect a weights array of complex exponentials (sometimes called twiddle factors) to exist prior to calling the function. Once created, this FFT weights array can be used over and over by the same Fourier function and can be shared by several Fourier functions.

FFT weights arrays are created by calling vDSP_create_fftsetup (single-precision) or vDSP_create_fftsetupD (double-precision). Before calling a function that processes in the frequency domain, you must call one of these functions. The caller specifies the required array size to the setup function. The setup functions:

• Create a data structure to hold the array.

• Build the array.

• Return a pointer to the data structure (or NULL if storage for this structure could not be allocated).

The pointer to the data structure is then passed as an argument to subsequent Fourier functions. You must check the value of the pointer before supplying the pointer to a frequency domain function. A returned value of zero is the only explicit notification that the array allocation failed.

Argument log2n to these functions is the base-2 logarithm of n, where n is the number of complex unit circle divisions the array represents, and thus specifies the largest number of elements that can be processed by a subsequent Fourier function. Argument log2n must equal or exceed argument log2n supplied to any functions using the weights array. Functions automatically adjust their strides through the array when the table has more resolution, or larger n, than required.

Thus, a single call to vDSP_create_fftsetup or vDSP_create_fftsetupD can serve a series of calls to FFT functions as long as n is large enough for each function called and as long as the weights array is preserved.

For example, if you need to call vDSP_fft_zrip for 256 data points, you must call vDSP_create_fftsetup with logn of at least 8. If you also need to call vDSP_fft_zrip for a set of 2048 data points, the value of logn must be at least 11. So you could generate this table once for both sets as follows:

FFTsetup setup; vDSP_create_fftsetup( 11, 0 ); /* supports up to 2048 (2**11) points  */

...

vDSP_fft_zrip( setup, small_set, 1, 8, FFT_FORWARD); /* 256 (2**8) points  */

...

vDSP_fft_zrip( setup, large_set, 1, 11, FFT_FORWARD); /* 2048 (2**11) points  */

Reusing a weights array for similarly sized FFTs is economical. However, using a weights array built for an FFT that processes a large number of elements can degrade performance for an FFT that processes a smaller number of elements.

For a more complete example of vDSP_create_fftsetup and FFT functions, see vDSP Examples.

Data Packing for Real FFTs

The discrete Fourier transform functions in the vDSP API provide a unique case in data formatting to conserve memory. Real-to-complex discrete Fourier transforms write their output data in special packed formats so that the complex output requires no more memory than the real input.

Packing for One-Dimensional Arrays

Due to its inherent symmetry, the forward Fourier transform of n floating-point inputs from the time domain to the frequency domain produces n/2 + 1 unique complex outputs. Because data point n/2 - i equals the complex conjugate of point n/2 + i, the first half of the frequency data by itself (up to and including point n/2 + 1) is sufficient to preserve the original information.

In addition, the FFT data packing takes advantage of the fact that the imaginary parts of the first complex output (element zero) and of the last complex output (element n/2, the Nyquist frequency), are always zero.

This makes it possible to store the real part of the last output where the imaginary part of the first output would have been. The zero-valued imaginary parts of the first and last outputs are implied.

For example, consider this eight-point real vector as it exists prior to being transformed to the frequency domain (Figure 2-1).

Transforming these eight real points to the frequency domain results in five complex values (Figure 2-2).

The five complex values are packed in the output vector shown in Figure 2-3

Notice that the real part of output four is stored where the imaginary part of output zero would have been stored with no packing. The imaginary parts of the first output and the last output, always zero, are implied.

Packing for Two-Dimensional Arrays

Two-dimensional real-to-complex DFTs are packed in a similar way. When processing two-dimensional real inputs, the DFTs first transform each row and then transform down columns. As each row is transformed, it is packed as previously described for one-dimensional transforms.

For example, to transform an 8-by-8 real image to the frequency domain, first each row is transformed and packed as shown in Figure 2-4

Notice that each row is packed like the one-dimensional transform in the earlier example. Packing yields 16 real values and 24 complex values in a total of 64 sizeof(float) memory locations, the size of the original matrix.

After each row is processed and packed, the DFT then transforms each column. Since the first pass leaves real values in columns zero and one, these two columns are transformed from real to complex using the same packing method vertically that was used horizontally on each row. Thus, transforming these 16 real values yields six complex points and four real points, occupying the 16 memory locations in columns zero and one.

The remaining three columns of complex values are transformed complex-to-complex, yielding 24 complex points as output, and requiring 48 memory locations to store. This brings the total storage requirement to 64 locations, the size of the original matrix. The completed 8-by-8 transform's format is shown in Figure 2-5

Scaling Fourier Transforms

To provide the best possible execution speeds, the vDSP library's functions don't always adhere strictly to textbook formulas for Fourier transforms, and must be scaled accordingly. The following sections specify the scaling for each type of Fourier transform implemented by the vDSP Library. The scaling factors are also stated explicitly in the formulas that accompany the function definitions in the reference chapter.

The scaling factors are summarized in Figure 2-6 which shows the implemented values in terms of the mathematical values. The sections that follow describe the scaling factors individually.

Real Fourier Transforms

Figure 2-7 shows the mathematical formula for the one-dimensional forward Fourier transform.

The values of the Fourier coefficients returned by the real forward transform as implemented (RF_imp) are equal to the mathematical values (RF_math) times 2:

Figure 2-8 shows the mathematical formula for the inverse Fourier transform.

The values of the Fourier coefficients returned by the real inverse transform as implemented (RI_imp) are equal to the mathematical values (RI_math) times N:

Complex Fourier Transforms

Figure 2-7 shows the mathematical formula for the one-dimensional forward Fourier transform.

The values of the Fourier coefficients returned by the complex forward transform as implemented (CF_imp) are equal to the mathematical values (CF_math):

The values of the Fourier coefficients returned by the complex inverse transform as implemented (CI_imp) are equal to the mathematical values (CI_math) times N:

Real 2-Dimensional Fourier Transforms

Figure 2-9 shows the mathematical formula for the 2-dimensional forward Fourier transform.

The values of the Fourier coefficients returned by the 2-dimensional real forward transform as implemented (RF2D_imp) are equal to the mathematical values (R2DF_math) times 2:

Figure 2-10 shows the mathematical formula for the 2-dimensional inverse Fourier transform.

The values of the Fourier coefficients returned by the 2-dimensional real inverse transform as implemented (RF2D_imp) are equal to the mathematical values (R2DF_math) times MN:

Complex 2-Dimensional Fourier Transforms

Figure 2-9 shows the mathematical formula for the 2-dimensional forward Fourier transform.

The values of the Fourier coefficients returned by the 2-dimensional complex forward transform as implemented (CF2D_imp) are equal to the mathematical values (CF2D_math):

Figure 2-10 shows the mathematical formula for the 2-dimensional inverse Fourier transform.

The values of the Fourier coefficients returned by the 2-dimensional complex inverse transform as implemented (CI2D_imp) are equal to the mathematical values (CI2D_math) times MN:

Scaling Example

Using complex two-dimensional transforms as an example, when transforming data from the time domain to the frequency domain with vDSP_fft_zop, no scaling factor is introduced. Transforming the image from the frequency domain to the time domain with vDSP_fft_zop, however, introduces a factor of N, where N is the vector length.

In the following example, a vector is transformed to the frequency domain, an inverse transform reconstructs it in the time domain, and then it is scaled by 1/N to recover the original input (perhaps altered by machine rounding errors).

Listing 2-1  Scaling example

    ~

    ~

scale = 1.0/((float) FFT_LENGTH) ;

    ~

    ~

/* to the frequency domain */

vDSP_fft_zop(   setup,

           &A,

           1,

           &B,

           1,

           FFT_LENGTH_LOG2,

           FFT_FORWARD,

           ) ;

/* back to the time domain */

vDSP_fft_zop(   setup,

           &B,

           1,

           &A,

           1,

           FFT_LENGTH_LOG2,

           FFT_INVERSE,

           ) ;

/*scale the result */

vDSP_vsmul(     A.realp,

           1,

           &scale,

           A.realp,

           1,

           FFT_LENGTH/2.0f

           ) ;

vDSP_vsmul(     A.imagp,

           1,

           &scale,

           A.imagp,

           1,

           FFT_LENGTH/2.0f

           ) ;

    ~

    ~