FFT/Fourier Transforms QuickStart Sample (F#)

Illustrates how to compute the forward and inverse Fourier transform of a real or complex signal using classes in the Extreme.Mathematics.SignalProcessing namespace in F#.

C# code Visual Basic code IronPython code Back to QuickStart Samples

// Illustrates the use of the FftProvider and Fft classes for computing 
// the Fourier transform of real and complex signals.

#light

open System

// We'll need real and complex vectors...
open Extreme.Mathematics
// The FFT classes reside in the Extreme.Mathematics.SignalProcessing
// namespace.
open Extreme.Mathematics.SignalProcessing

// This QuickStart sample shows how to compute the Fouier
// transform of real and complex signals.

// Some vectors to play with:
let r1 = Vector.Create(1000, fun i -> 1.0 / (float (1+i)))
let c1 = ComplexVector.Create(1000, fun i -> new DoubleComplex(Math.Sin(0.03 * float i), Math.Cos(0.07 * float i)))

let r2 = Vector.Create(1.0, 2.0, 3.0, 4.0)
let c2 = 
    ComplexVector.Create(
        new DoubleComplex(1.0, 2.0), 
        new DoubleComplex(3.0, 4.0), 
        new DoubleComplex(5.0, 6.0), 
        new DoubleComplex(7.0, 8.0))

//
// One-time FFT's
//

// The Vector and ComplexVector classes have static methods to compute FFT's:
let c3 = Vector.FourierTransform(r2)
let r3 = Vector.InverseFourierTransform(c3)
printfn "fft(r2) = %s" (c3.ToString("F3"))
printfn "ifft(fft(r2)) = %s" (r3.ToString("F3"))
// The ComplexConjugateSignalVector type represents a complex vector
// that is the Fourier transform of a real signal. 
// It enforces certain symmetry properties:
printfn "c3.[i] == conj(c3.[N-i]): %A == conj(%A)" c3.[1] c3.[3]

//
// FFT Providers
//

// FFT's require a fair bit of pre-computation. Using the FftProvider class,
// you can get an Fft object that caches these computations.

// Here, we create an FFT implementation for a real signal:
let realFft = FftProvider.ManagedProvider.Create1DRealFft(r1.Length)
// For a complex to complex transform:
let complexFft = FftProvider.ManagedProvider.Create1DComplexFft(c1.Length)

// You can set the scale factor for the forward transform.
// The default is 1/N.
realFft.ForwardScaleFactor <- 1.0 / Math.Sqrt(float c1.Length)
// and the backward transform, with default 1:
realFft.BackwardScaleFactor <- realFft.ForwardScaleFactor

// The ForwardTransform method performs a forward transform:
let c4 = realFft.ForwardTransform(r1)
printfn "First 5 terms of fft(r1):"
for i in 0..4 do
    printfn "   %d: %A" i c4.[i]
let c5 = complexFft.ForwardTransform(c1)
printfn "First 5 terms of fft(c1):"
for i in 0..4 do
    printfn "   %d: %A" i c5.[i]

// ForwardTransform has many overloads for real to complex and
// complex to complex transforms.

// A one-sided transform returns only the first half of the FFT of
// a real signal. The rest can be deduced from the symmetry properties.
// Here's how to compute a one-sided FFT:
let c6 = ComplexVector.Create(r1.Length / 2 + 1)
realFft.ForwardTransform(r1, c6, RealFftFormat.OneSided)

// The BackwardTransform method has a similar set of overloads:
let r4 = Vector.Create(r1.Length)
realFft.BackwardTransform(c6, r4, RealFftFormat.OneSided)

// As the last step, we need to dispose the two FFT implementations.
realFft.Dispose()
complexFft.Dispose()

//
// 2D transforms
//

// 2D transforms are handled in a completely analogous way.
let m = Matrix.Create(36, 56, fun i -> fun j -> exp(-0.1 * float i) * Math.Sin(0.01 * (float (i * i + j * j - i * j))))
let mFft = ComplexMatrix.Create(m.RowCount, m.ColumnCount)

let fft2 = FftProvider.CurrentProvider.Create2DRealFft(m.RowCount, m.ColumnCount)
fft2.ForwardTransform(m, mFft)

printfn "First few terms of fft(m):"
for i in 0..3 do
    for j in 0..3 do 
        if j > 0 then 
            printf ", "
        printf "%s" (mFft.[i, j].ToString("F4"))
    printfn ""

// and the backward transform:
fft2.BackwardTransform(mFft, m)
fft2.Dispose()

// Dispose is called automatically.

printf "Press Enter key to exit..."
Console.ReadLine() |> ignore