Fast Fourier Transform

The fft package is a wrapper of the C language version of PocketFFT library designed to support FFT of real to complex and complex to real (arrays).

The result of a real-to-complex transform, because of mathematical symmetry of the result, is stored in the original input array rather than 2x the space.

The output is the two real values bracketing the complex pairs of conjugate negative frequencies: r0 r1 i1 r2 i2 r3 i3 ... rx

where r0 + i0 is the first complex result, r1 - i1 is the second, and so on until rx + i0 (where x is n/2) is the last. (Note the minus signs.)

The positive frequencies are the same as the negative frequencies in reverse order. See the reference for FFTW for further examples of embeddings.

Interleaved complex transforms

The complex FFT API supports both f32 and f64 interleaved buffers in the form [real0, imag0, real1, imag1, ...]. Create a plan with create_complex_plan_f32(length) or create_complex_plan_f64(length), then call the matching forward_complex_f32/backward_complex_f32 or forward_complex_f64/backward_complex_f64 functions. They transform in place, return 0 on success, and return -1 when the plan or buffer length does not match. The inverse is unnormalized; divide each component by the original transform length for the normalized result. Destroy each plan exactly once with destroy_plan.

The helpers require exactly twice the plan length in values of the matching precision and validate that size before calling the native backend.

import vsl.fft

mut samples := [f64(1), 0, 0, 0] // two complex values
plan := fft.create_complex_plan_f64(2)!
defer {
    fft.destroy_plan(plan)
}
fft.forward_complex_f64(plan, mut samples)
// Divide inverse results by 2 to normalize them.
fft.backward_complex_f64(plan, mut samples)

The single-precision API has the same plan and buffer contract:

import vsl.fft

mut samples := [f32(1), 0, 0, 0]
plan := fft.create_complex_plan_f32(2)!
defer {
    fft.destroy_plan(plan)
}
fft.forward_complex_f32(plan, mut samples)
fft.backward_complex_f32(plan, mut samples)

fn backward_complex_f32(plan Fftplan, mut data []f32) int

backward_complex_f32 applies the unnormalized inverse transform to interleaved complex f32 values in place. Divide every result by the transform length to obtain the normalized inverse.

fn backward_complex_f64(plan Fftplan, mut data []f64) int

backward_complex_f64 applies the unnormalized inverse transform to interleaved complex f64 values in place. The data slice must contain exactly 2 * plan_length values. Divide every result by the transform length to obtain the normalized inverse.

fn backward_fft[T](r Fftplan, mut v T) int

backward_fft computes the backwards Fourier transform defined by the plan r.

fn create_complex_plan_f32(length int) !Fftplan

create_complex_plan_f32 creates a reusable single-precision plan for interleaved complex data stored as [real0, imag0, real1, imag1, ...].

fn create_complex_plan_f64(length int) !Fftplan

create_complex_plan_f64 creates a reusable double-precision plan for interleaved complex data stored as [real0, imag0, real1, imag1, ...].

fn create_plan[T](x T) ?Fftplan

create_plan returns a plan to compute a Fourier transform of the given array. A plan is reusable for any array of exactly this size and type. The array may be []f32, []f64, or []complx_f32 or []complex_f64.

fn destroy_plan(plan Fftplan)

destroy_plan releases the native memory owned by a plan returned from create_plan. Each plan must be destroyed exactly once after its last use.

fn forward_complex_f32(plan Fftplan, mut data []f32) int

forward_complex_f32 transforms interleaved complex f32 values in place. The data slice must contain exactly 2 * plan_length values.

fn forward_complex_f64(plan Fftplan, mut data []f64) int

forward_complex_f64 transforms interleaved complex f64 values in place. The data slice must contain exactly 2 * plan_length values.

fn forward_fft[T](p Fftplan, mut v T) int

forward_fft computes a Fourier transform defined by the plan p. The input is []f32 or []f64. The output result (r) is returned in-place, and is complex and real numbers mixed as follows: r[0] (im[0] is assumed and is 0) r[1] im[1] ... r[n/2] im[n/2] conjugate negative frequences r[n] is the magnitude at zero frequency (im[n] is assumed and is 0) Positive frequencies are excluded as these are symmetric to the negative frequencies for real inputs. (Consult any signal processing text for details.)

To generate a complete result, copy the results to a larger complex array using r[0] and i0 r[1] and i-r[2] r[3] and i*-r[4] ... r[n] and i0 at the zero frequency r[n-2] and ir[n-1] the complex conjugate negative frequencies r[n-4] and ir[n-3] in reverse order ... r[1] and ir[2]

Note: these codes allocate and free memory of the same size as the input.