copy code from original project

benchmarks/Benchmarks.elm

@@ -0,0 +1,91 @@
+module BenchmarksMain exposing (..)
+
+import Array
+import Benchmark exposing (Benchmark, benchmark, describe, scale)
+import Benchmark.Runner exposing (program)
+import FFT.FFT as FFT exposing (real)
+import Random
+
+
+main =
+    program fftbenchmarks
+
+
+fftbenchmarks : Benchmark
+fftbenchmarks =
+    let
+        input128 =
+            List.range 0 127
+                |> List.map (toFloat >> (/) 128 >> (*) pi >> cos >> real)
+
+        input128NoiseFloat =
+            Random.step
+                (Random.list 128 (Random.float -1 1))
+                (Random.initialSeed 154324632)
+                |> Tuple.first
+
+        input128NoiseComplex =
+            input128NoiseFloat
+                |> List.map real
+
+        input256NoiseFloat =
+            input128NoiseFloat ++ input128NoiseFloat
+
+        input512NoiseFloat =
+            input256NoiseFloat ++ input256NoiseFloat
+
+        input256 =
+            List.range 0 255
+                |> List.map (toFloat >> (/) 128 >> (*) pi >> cos >> real)
+
+        input512 =
+            List.range 0 511
+                |> List.map (toFloat >> (/) 128 >> (*) pi >> cos >> real)
+
+        input1024 =
+            List.append input512 input512
+
+        inputArray128 =
+            Array.fromList input128
+
+        inputArray256 =
+            Array.fromList input256
+
+        inputArray512 =
+            Array.fromList input512
+
+        inputArray1024 =
+            Array.append inputArray512 inputArray512
+
+        makeInput index =
+            ( (2 ^ toFloat index |> round |> String.fromInt) ++ " samples"
+            , List.range 0 (2 ^ toFloat index |> round)
+                |> List.map (toFloat >> (/) 128 >> (*) pi >> cos >> real)
+            )
+    in
+    describe "FFT"
+        [ --[ scale "list based" <|
+          --    (List.range 5 8
+          --        |> List.map makeInput
+          --        |> List.map (Tuple.mapSecond (\input -> \() -> List.length input))
+          --     --|> List.map (Tuple.mapSecond (\input -> \() -> FFT.fft input))
+          --    )
+          --,
+          describe "convolve"
+            [ Benchmark.benchmark "128 samples"
+                (\() ->
+                    FFT.convolve input128NoiseFloat input128NoiseFloat
+                )
+            , Benchmark.benchmark "256 samples"
+                (\() ->
+                    FFT.convolve input256NoiseFloat input256NoiseFloat
+                )
+            ]
+
+        --,
+        --Benchmark.compare "128 samples List vs Array"
+        --  "List"
+        --  (\_ -> FFT.fft input128)
+        --  "Array"
+        --  (\_ -> FFT.fftArray inputArray128)
+        ]

elm.json

@@ -9,9 +9,11 @@
     ],
     "elm-version": "0.19.0 <= v < 0.20.0",
     "dependencies": {
-        "elm/core": "1.0.2 <= v < 2.0.0"
+        "elm/core": "1.0.2 <= v < 2.0.0",
+        "elm-community/list-extra": "8.7.0 <= v < 9.0.0",
+        "elm-explorations/benchmark": "1.0.2 <= v < 2.0.0"
     },
     "test-dependencies": {
-        "elm-explorations/test": "1.0.0 <= v < 2.0.0"
+        "elm-explorations/test": "2.0.0 <= v < 3.0.0"
     }
-}
+}

src/Complex.elm

@@ -0,0 +1,113 @@
+module Complex exposing (..)
+
+{-| Simple representation of a complex number
+-}
+
+
+type alias Complex =
+    { re : Float
+    , im : Float
+    }
+
+
+{-| Construct a `Complex` number with imaginary value only
+
+    imaginary 1 == { re = 0, im = 1 }
+
+-}
+imaginary : Float -> Complex
+imaginary im =
+    { re = 0, im = im }
+
+
+{-| Construct a `Complex` number with real value only
+
+    real 1 == { re = 1, im = 0 }
+
+-}
+real : Float -> Complex
+real re =
+    { re = re, im = 0 }
+
+
+{-| The complex zero
+-}
+zero : Complex
+zero =
+    { re = 0, im = 0 }
+
+
+{-| Check equality of two complex numbers
+-}
+equal : Complex -> Complex -> Bool
+equal a b =
+    a.re == b.re && a.im == b.im
+
+
+{-| Add two complex numbers
+
+    add (real 1) (imaginary 2) == { re = 1, im = 2 }
+
+-}
+add : Complex -> Complex -> Complex
+add a b =
+    { re = a.re + b.re, im = a.im + b.im }
+
+
+{-| Subtract two complex numbers
+
+    subtract { re = 1, im = 2 } (real 2) == { re = -1, im = 2 }
+
+-}
+subtract : Complex -> Complex -> Complex
+subtract a b =
+    { re = a.re - b.re, im = a.im - b.im }
+
+
+{-| Multiply two complex numbers
+
+    multiply (imaginary 1) (real 2) == { re = 0, im = 2 }
+
+-}
+multiply : Complex -> Complex -> Complex
+multiply a b =
+    { re = (a.re * b.re) - (a.im * b.im)
+    , im = (a.re * b.im) + (a.im * b.re)
+    }
+
+
+{-| Divide a complex number by a real number.
+Notice that the divisor is given first.
+
+    real 2 |> divideByReal 2 == { re = 0, im = 2 }
+
+-}
+divideByReal : Float -> Complex -> Complex
+divideByReal divisor z =
+    { re = z.re / divisor, im = z.im / divisor }
+
+
+{-| Given complex number z, computes e^z
+-}
+exp : Complex -> Complex
+exp { re, im } =
+    { re = (Basics.e ^ re) * cos im
+    , im = (e ^ re) * sin im
+    }
+
+
+{-| The absolute value of complex number, which is also its distance from zero
+-}
+abs : Complex -> Float
+abs { re, im } =
+    sqrt ((re ^ 2) + (im ^ 2))
+
+
+{-| computes the conjugate of a complex number
+
+    conjugate { re = 1, im = 1 } == { re = 1, im = -1 }
+
+-}
+conjugate : Complex -> Complex
+conjugate z =
+    { re = z.re, im = -z.im }

src/Convolution.elm

@@ -0,0 +1,32 @@
+module Convolution exposing (..)
+
+import Complex exposing (multiply)
+import FFT.FFT exposing (inverseFft, realFft)
+import List.Extra
+
+
+convolve : List Float -> List Float -> List Float
+convolve xs ys =
+    if List.length xs >= 128 || List.length ys >= 128 then
+        let
+            fxs =
+                realFft xs
+
+            fys =
+                realFft ys
+        in
+        List.map2 multiply fxs fys
+            |> inverseFft
+            |> List.map .re
+
+    else
+        convolveDirectly xs ys
+
+
+convolveDirectly : List Float -> List Float -> List Float
+convolveDirectly xs ys =
+    List.Extra.initialize (List.length xs)
+        (\n ->
+            List.map2 (*) xs (List.drop n ys)
+                |> List.sum
+        )

src/FFT.elm

@@ -0,0 +1,107 @@
+module FFT.FFT exposing (fft, inverseFft, realFft, realInverseFft, zeroPadTo, zeroPadToNextPowerOf2)
+
+import Complex exposing (..)
+
+
+zeroPadTo : Int -> List Complex -> List Complex
+zeroPadTo n xs =
+    if List.length xs >= n then
+        xs
+
+    else
+        xs ++ List.repeat (n - List.length xs) zero
+
+
+zeroPadToNextPowerOf2 : List Complex -> List Complex
+zeroPadToNextPowerOf2 zs =
+    let
+        length =
+            List.length zs
+
+        nextPowerOf2 : Int
+        nextPowerOf2 =
+            2 ^ (ceiling <| logBase 2 (toFloat length))
+    in
+    zeroPadTo nextPowerOf2 zs
+
+
+splitEvensAndOdds : List a -> ( List a, List a )
+splitEvensAndOdds list =
+    case list of
+        [] ->
+            ( [], [] )
+
+        x :: [] ->
+            ( [ x ], [] )
+
+        x :: y :: xs ->
+            let
+                ( xs_, ys_ ) =
+                    splitEvensAndOdds xs
+            in
+            ( x :: xs_, y :: ys_ )
+
+
+fft : List Complex -> List Complex
+fft array =
+    let
+        fftInner xs =
+            case List.length xs of
+                0 ->
+                    []
+
+                1 ->
+                    xs
+
+                n ->
+                    let
+                        ( evens, odds ) =
+                            splitEvensAndOdds xs
+
+                        ys =
+                            fft evens
+
+                        zs =
+                            fft odds
+
+                        ts =
+                            List.indexedMap
+                                (\i z ->
+                                    multiply (coefficient i) z
+                                )
+                                zs
+
+                        coefficient k =
+                            exp (imaginary <| -2 * pi * (toFloat k / toFloat n))
+                    in
+                    List.map2 add ys ts
+                        ++ List.map2 subtract ys ts
+    in
+    fftInner (zeroPadToNextPowerOf2 array)
+
+
+inverseFft : List Complex -> List Complex
+inverseFft list =
+    let
+        paddedInput =
+            zeroPadToNextPowerOf2 list
+    in
+    paddedInput
+        |> List.map conjugate
+        |> fft
+        |> List.map conjugate
+        |> List.map (divideByReal (List.length paddedInput |> toFloat))
+
+
+realFft : List Float -> List Complex
+realFft floats =
+    floats
+        |> List.map real
+        |> fft
+
+
+realInverseFft : List Float -> List Complex
+realInverseFft floats =
+    floats
+        |> List.map real
+        |> inverseFft