3 people like it.

Optimized Fast Fourier Transform (FFT)

Fast Fourier Transform algorithm inspired by the snippet http://fssnip.net/dC/title/fast-Fourier-transforms-FFT-. The original snippet was a beautiful and idiomatic implementation based on the definition of FFT. This version sacrifices idioms and elegance for performance. For data sets with 1024 samples this version seems to be about 100x faster and easier on the GC.

  1: 
  2: 
  3: 
  4: 
  5: 
  6: 
  7: 
  8: 
  9: 
 10: 
 11: 
 12: 
 13: 
 14: 
 15: 
 16: 
 17: 
 18: 
 19: 
 20: 
 21: 
 22: 
 23: 
 24: 
 25: 
 26: 
 27: 
 28: 
 29: 
 30: 
 31: 
 32: 
 33: 
 34: 
 35: 
 36: 
 37: 
 38: 
 39: 
 40: 
 41: 
 42: 
 43: 
 44: 
 45: 
 46: 
 47: 
 48: 
 49: 
 50: 
 51: 
 52: 
 53: 
 54: 
 55: 
 56: 
 57: 
 58: 
 59: 
 60: 
 61: 
 62: 
 63: 
 64: 
 65: 
 66: 
 67: 
 68: 
 69: 
 70: 
 71: 
 72: 
 73: 
 74: 
 75: 
 76: 
 77: 
 78: 
 79: 
 80: 
 81: 
 82: 
 83: 
 84: 
 85: 
 86: 
 87: 
 88: 
 89: 
 90: 
 91: 
 92: 
 93: 
 94: 
 95: 
 96: 
 97: 
 98: 
 99: 
100: 
101: 

// For examples and tests see: https://gist.github.com/mrange/da57f972b3dfdfb44f28fd340841586c
//  Inspired by: http://fssnip.net/dC/title/fast-Fourier-transforms-FFT-
//  and: https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
//  Sacrifices idioms for performance
module FourierTransform =
  open System
  open System.Numerics

  let maxSize       = 4096
  let pi            = Math.PI
  let tau           = 2.*pi

  module internal Details =
    let isPowerOf2 n  = (n &&& n - 1) = 0
    let ilog2 n       =
      if n < 2              then failwith "n must be greater than 1"
      if not (isPowerOf2 n) then failwith "n must be a power of 2"
      let rec loop n c s =
        let t = 1 <<< c
        if t = n then
          c
        elif t > n then
          loop n (c - s) (s >>> 1)
        else
          loop n (c + s) (s >>> 1)
      loop n 16 8
    let twiddle a     = Complex.FromPolarCoordinates(1., a)

    let twiddles      =
      let unfolder c =
        if c < 2*maxSize then
          let vs = Array.init (c / 2) (fun i -> twiddle (-tau * float i / float c))
          Some (vs, c*2)
        else
          None
      Array.unfold unfolder 1

    let rec loop n2 ln s c f t =
      if c > 2 then
        let c2 = c >>> 1
        let struct (t, f) = loop n2 (ln - 1) (s <<< 1) c2 f t

        let twiddles = twiddles.[ln]

        if s > 1 then
          for j = 0 to c2 - 1 do
            let w   = twiddles.[j]
            let off = s*j
            let off2= off <<< 1;
            for i = 0 to s - 1 do
              let e = Array.get f (i + off2 + 0)
              let o = Array.get f (i + off2 + s)
              let a = w*o
              Array.set t (i + off + 0)   (e + a)
              Array.set t (i + off + n2)  (e - a)
        else
          for j = 0 to c2 - 1 do
            let w = twiddles.[j]
            let e = Array.get f (2*j + 0)
            let o = Array.get f (2*j + s)
            let a = w*o
            Array.set t (j + 0)   (e + a)
            Array.set t (j + n2)  (e - a)

        struct (f, t)
      elif c = 2 then
        for i = 0 to s - 1 do
          let e = Array.get f (i + 0)
          let o = Array.get f (i + s)
          let a = o
          Array.set t (i + 0)   (e + a)
          Array.set t (i + n2)  (e - a)

        struct (f, t)
      else
        struct (t, f)

  open Details
  
  let dft (vs : Complex []) =
    let l   = vs.Length
    let am  = tau / float l
    let rec loop s j i =
      if j < l then
        let v = vs.[j]
        let n = v*twiddle (-float i * float j * am)
        loop (s + n) (j + 1) i
      else 
        s
    Array.init l (loop Complex.Zero 0)

  let fft (vs : Complex []) : Complex [] =
    let n   = vs.Length
    let ln  = ilog2 n

    let vs0 = Array.copy vs
    let vs1 = Array.zeroCreate n

    let struct (_, t) = Details.loop (n >>> 1) ln 1 n vs0 vs1

    t
namespace System
namespace System.Numerics
val maxSize : int

Full name: Script.FourierTransform.maxSize
val pi : float

Full name: Script.FourierTransform.pi
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
field Math.PI = 3.14159265359
val tau : float

Full name: Script.FourierTransform.tau
module Details

from Script.FourierTransform
val internal isPowerOf2 : n:int -> bool

Full name: Script.FourierTransform.Details.isPowerOf2
val n : int
val internal ilog2 : n:int -> int32

Full name: Script.FourierTransform.Details.ilog2
val failwith : message:string -> 'T

Full name: Microsoft.FSharp.Core.Operators.failwith
val not : value:bool -> bool

Full name: Microsoft.FSharp.Core.Operators.not
val loop : (int -> int32 -> int32 -> int32)
val c : int32
val s : int32
val t : int
val internal twiddle : a:float -> Complex

Full name: Script.FourierTransform.Details.twiddle
val a : float
Multiple items
type Complex =
  struct
    new : real:float * imaginary:float -> Complex
    member Equals : obj:obj -> bool + 1 overload
    member GetHashCode : unit -> int
    member Imaginary : float
    member Magnitude : float
    member Phase : float
    member Real : float
    member ToString : unit -> string + 3 overloads
    static val Zero : Complex
    static val One : Complex
    ...
  end

Full name: System.Numerics.Complex

--------------------
Complex()
Complex(real: float, imaginary: float) : unit
Complex.FromPolarCoordinates(magnitude: float, phase: float) : Complex
val internal twiddles : Complex [] []

Full name: Script.FourierTransform.Details.twiddles
val unfolder : (int -> (Complex [] * int) option)
val c : int
val vs : Complex []
type Array =
  member Clone : unit -> obj
  member CopyTo : array:Array * index:int -> unit + 1 overload
  member GetEnumerator : unit -> IEnumerator
  member GetLength : dimension:int -> int
  member GetLongLength : dimension:int -> int64
  member GetLowerBound : dimension:int -> int
  member GetUpperBound : dimension:int -> int
  member GetValue : [<ParamArray>] indices:int[] -> obj + 7 overloads
  member Initialize : unit -> unit
  member IsFixedSize : bool
  ...

Full name: System.Array
val init : count:int -> initializer:(int -> 'T) -> 'T []

Full name: Microsoft.FSharp.Collections.Array.init
val i : int
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
union case Option.Some: Value: 'T -> Option<'T>
union case Option.None: Option<'T>
val unfold : generator:('State -> ('T * 'State) option) -> state:'State -> 'T []

Full name: Microsoft.FSharp.Collections.Array.unfold
val internal loop : n2:int -> ln:int -> s:int -> c:int -> f:Complex [] -> t:Complex [] -> unit

Full name: Script.FourierTransform.Details.loop
val n2 : int
val ln : int
val s : int
val f : Complex []
val t : Complex []
val c2 : int
val twiddles : Complex []
val j : int
val w : Complex
val off : int
val off2 : int
val e : Complex
val get : array:'T [] -> index:int -> 'T

Full name: Microsoft.FSharp.Collections.Array.get
val o : Complex
val a : Complex
val set : array:'T [] -> index:int -> value:'T -> unit

Full name: Microsoft.FSharp.Collections.Array.set
module Array

from Microsoft.FSharp.Collections
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = System.Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
val copy : array:'T [] -> 'T []

Full name: Microsoft.FSharp.Collections.Array.copy
val zeroCreate : count:int -> 'T []

Full name: Microsoft.FSharp.Collections.Array.zeroCreate
Raw view Test code New version

More information

Link:http://fssnip.net/7Tn
Posted:6 years ago
Author:Mårten Rånge
Tags: algorithm , performance , fft

This web site is created using F# and Suave web server. It is hosted on Azure and the source code is on GitHub. Contributions are welcome!

The first version of fssnip.net has been created by @tomaspetricek back in 2010. This web site is a new, open-source and contribution-friendly version.