first

README.md

@@ -0,0 +1,4 @@
+# vector-fft
+
+Native FFT and IFFT for vector
+

Setup.hs

@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain

src/Data/Vector/FFT.hs

@@ -0,0 +1,150 @@
+{-# language BangPatterns        #-}
+{-# language LambdaCase          #-}
+
+module Data.Vector.FFT where
+
+import Control.Monad (when)
+import Control.Monad.Primitive (PrimMonad(..))
+
+import Control.Monad.ST (runST)
+import Data.Bits (shiftR,shiftL,(.&.),(.|.))
+import Data.Bool (Bool,otherwise)
+import Data.Complex (Complex(..),conjugate)
+
+import Data.Vector.Unboxed as V (Vector, Unbox(..), map, length, unsafeFreeze)
+import qualified Data.Vector.Unboxed.Mutable as VM (read, write, copy, new)
+
+
+import Prelude hiding (read)
+
+{-# RULES
+"fft/ifft" forall x. fft (ifft x) = x
+"ifft/fft" forall x. ifft (fft x) = x
+  #-}
+
+-- | Radix-2 decimation-in-time fast Fourier Transform.
+--   The given array must have a length that is a power of two.
+
+{-# inlinable [1] fft #-}
+fft arr = if arrOK arr
+  then runST $ do {
+      marr <- copyWhole arr
+    ; mfft marr
+    ; V.unsafeFreeze marr
+  }
+  else Prelude.error "Data.Vector.FFT.fft: bad array length"
+
+-- | Inverse fast Fourier transform.
+
+{-# inlinable [1] ifft #-}
+ifft arr = if arrOK arr
+  then
+    let lenComplex = intToComplexDouble (V.length arr)
+    in cmap ((/ lenComplex) . conjugate) . fft . cmap conjugate $ arr
+  else error "Data.Vector.FFT.ifft: bad vector length"
+
+
+{-# inline copyWhole #-}
+copyWhole arr = do
+  let len = V.length arr
+  marr <- VM.new len
+  VM.copy marr arr
+  pure marr
+
+
+{-# inline arrOK #-}
+arrOK arr =
+  let n = V.length arr
+  in (1 `shiftL` log2 n) == n
+
+-- | Radix-2 decimation-in-time fast Fourier Transform.
+--   The given array must have a length that is a power of two,
+--   though this property is not checked.
+
+mfft mut = do {
+    let len = V.length mut
+  ; let bitReverse !i !j = do {
+          ; if i == len - 1
+              then stage 0 1
+              else do {
+                  when (i < j) $ swap mut i j
+                ; let inner k l = if k <= l
+                        then inner (k `shiftR` 1) (l - k)
+                        else bitReverse (i + 1) (l + k)
+                ; inner (len `shiftR` 1) j
+              }
+        }
+        stage l l1 = if l == (log2 len)
+          then pure ()
+          else do {
+              let !l2 = l1 `shiftL` 1
+                  !e = (negate twoPi) / (intToDouble l2)
+                  flight j !a = if j == l1
+                    then stage (l + 1) l2
+                    else do {
+                        let butterfly i = if i >= len
+                              then flight (j + 1) (a + e)
+                              else do {
+                                  let i1 = i + l1
+                                ; xi1 :+ yi1 <- VM.read mut i1
+                                ; let !c = cos a
+                                      !s = sin a
+                                      d = (c*xi1 - s*yi1) :+ (s*xi1 + c*yi1)
+                                ; ci <- VM.read mut i
+                                ; VM.write mut i1 (ci - d)
+                                ; VM.write mut i (ci + d)
+                                ; butterfly (i + l2)
+                              }
+                      ; butterfly j
+                    }
+            ; flight 0 0
+         }
+  ; bitReverse 0 0
+}
+
+-- wildcard cases should never happen. if they do, really bad things will happen.
+b,s :: Int -> Int
+b = \case { 0 -> 0x02; 1 -> 0x0c; 2 -> 0xf0; 3 -> 0xff00; 4 -> wordToInt 0xffff0000; 5 -> wordToInt 0xffffffff00000000; _ -> 0; }
+s = \case { 0 -> 1; 1 -> 2; 2 -> 4; 3 -> 8; 4 -> 16; 5 -> 32; _ -> 0; }
+{-# inline b #-}
+{-# inline s #-}
+
+log2 :: Int -> Int
+log2 v0 = if v0 <= 0
+  then error $ "Data.Vector.FFT: nonpositive input, got " ++ show v0
+  else go 5 0 v0
+  where
+    go !i !r !v
+      | i == -1 = r
+      | v .&. b i /= 0 =
+          let si = s i
+          in go (i - 1) (r .|. si) (v `shiftR` si)
+      | otherwise = go (i - 1) r v
+
+
+{-# inline swap #-}
+swap mut i j = do
+  atI <- VM.read mut i
+  atJ <- VM.read mut j
+  VM.write mut i atJ
+  VM.write mut j atI
+
+twoPi :: Double
+{-# inline twoPi #-}
+twoPi = 6.283185307179586
+
+intToDouble :: Int -> Double
+{-# inline intToDouble #-}
+intToDouble = fromIntegral
+
+wordToInt :: Word -> Int
+{-# inline wordToInt #-}
+wordToInt = fromIntegral
+
+intToComplexDouble :: Int -> Complex Double
+{-# inline intToComplexDouble #-}
+intToComplexDouble = fromIntegral
+
+
+{-# inline cmap #-}
+cmap = V.map

stack.yaml

@@ -0,0 +1,67 @@
+# This file was automatically generated by 'stack init'
+#
+# Some commonly used options have been documented as comments in this file.
+# For advanced use and comprehensive documentation of the format, please see:
+# https://docs.haskellstack.org/en/stable/yaml_configuration/
+
+# Resolver to choose a 'specific' stackage snapshot or a compiler version.
+# A snapshot resolver dictates the compiler version and the set of packages
+# to be used for project dependencies. For example:
+#
+# resolver: lts-3.5
+# resolver: nightly-2015-09-21
+# resolver: ghc-7.10.2
+#
+# The location of a snapshot can be provided as a file or url. Stack assumes
+# a snapshot provided as a file might change, whereas a url resource does not.
+#
+# resolver: ./custom-snapshot.yaml
+# resolver: https://example.com/snapshots/2018-01-01.yaml
+resolver:
+  url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/16/27.yaml
+
+# User packages to be built.
+# Various formats can be used as shown in the example below.
+#
+# packages:
+# - some-directory
+# - https://example.com/foo/bar/baz-0.0.2.tar.gz
+#   subdirs:
+#   - auto-update
+#   - wai
+packages:
+- .
+# Dependency packages to be pulled from upstream that are not in the resolver.
+# These entries can reference officially published versions as well as
+# forks / in-progress versions pinned to a git hash. For example:
+#
+# extra-deps:
+# - acme-missiles-0.3
+# - git: https://github.com/commercialhaskell/stack.git
+#   commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a
+#
+# extra-deps: []
+
+# Override default flag values for local packages and extra-deps
+# flags: {}
+
+# Extra package databases containing global packages
+# extra-package-dbs: []
+
+# Control whether we use the GHC we find on the path
+# system-ghc: true
+#
+# Require a specific version of stack, using version ranges
+# require-stack-version: -any # Default
+# require-stack-version: ">=2.5"
+#
+# Override the architecture used by stack, especially useful on Windows
+# arch: i386
+# arch: x86_64
+#
+# Extra directories used by stack for building
+# extra-include-dirs: [/path/to/dir]
+# extra-lib-dirs: [/path/to/dir]
+#
+# Allow a newer minor version of GHC than the snapshot specifies
+# compiler-check: newer-minor

test/Spec.hs

@@ -0,0 +1 @@
+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}

vector-fft.cabal

@@ -0,0 +1,47 @@
+name:                vector-fft
+version:             0.1.0.0
+synopsis:            Native FFT and IFFT for vector
+-- description:
+homepage:            https://github.com/ocramz/vector-fft
+license:             BSD3
+license-file:        LICENSE
+author:              chessai, Marco Zocca
+maintainer:          ocramz
+copyright:           2020 chessai, ocramz
+category:            Numeric
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+tested-with:         GHC == 7.10.2
+
+library
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  hs-source-dirs:      src
+  exposed-modules:     Data.Vector.FFT
+  build-depends:       base >= 4.7 && < 5
+                     , primitive
+                     , vector
+
+-- executable vector-fft
+--   default-language:    Haskell2010
+--   ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+--   hs-source-dirs:      app
+--   main-is:             Main.hs
+--   build-depends:       base
+--                      , vector-fft
+
+test-suite spec
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  build-depends:       base
+                     , vector-fft
+                     , hspec
+                     , QuickCheck
+
+source-repository head
+  type:     git
+  location: https://github.com/githubuser/vector-fft