contiguous-fft 0.1.0.1 → 0.2.0.0

raw patch · 2 files changed

+178/−148 lines, 2 filesdep +semiringsdep −prim-instancesdep −primitivedep ~basedep ~contiguousPVP ok

version bump matches the API change (PVP)

Dependencies added: semirings

Dependencies removed: prim-instances, primitive

Dependency ranges changed: base, contiguous

API changes (from Hackage documentation)

Files

• contiguous-fft.cabal +22/−2
• src/Data/Primitive/Contiguous/FFT.hs +156/−146

@@ -1,5 +1,5 @@ name:                contiguous-fft-version:             0.1.0.1+version:             0.2.0.0 synopsis:            dft of contiguous memory structures description:         DFT and iDFT on data structures implementing a common                      contiguous interface@@ -14,9 +14,29 @@ extra-source-files:  ChangeLog.md cabal-version:       >=1.10 +source-repository head+  type: git+  location: https://github.com/chessai/contiguous-fft.git+ library   exposed-modules:     Data.Primitive.Contiguous.FFT-  build-depends:       base >=4.9 && <5.0, contiguous >=0.2.0.0, prim-instances, primitive >= 0.6.4.0+  build-depends:+      base >=4.9 && <5+    , contiguous >=0.3+    , semirings >= 0.3   hs-source-dirs:      src   default-language:    Haskell2010+  ghc-options: -O2++-- test-suite spec+--   type: exitcode-stdio-1.0+--   hs-source-dirs: test+--   main-is: Spec.hs+--   build-depends:+--      base >= 4.9 && < 5+--    , primitive >= 0.6.4.0+--    , hedgehog >= 0.6+--    , contiguous-fft+--    , prim-instances+--    , math-functions

@@ -1,166 +1,176 @@-{-# LANGUAGE BangPatterns        #-}-{-# LANGUAGE NoImplicitPrelude   #-}-{-# LANGUAGE ScopedTypeVariables #-}+{-# language BangPatterns        #-}+{-# language LambdaCase          #-}+{-# language NoImplicitPrelude   #-}+{-# language ScopedTypeVariables #-} +-- | This module exposes functions for performing+--   Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT)+--   over 'Contiguous' data structures. module Data.Primitive.Contiguous.FFT-  ( dft-  , idft-  , dftMutable-  , overlapDFT+  ( fft+  , ifft   ) where  import qualified Prelude -import Data.Eq (Eq((==)))-import Data.Function (($))-import Control.Monad-import Data.Ord-import Control.Monad.ST-import Data.Complex hiding (cis)-import qualified Data.Complex as C-import Data.Primitive.Contiguous-import GHC.Num (Num(..))-import GHC.Float-import GHC.Real-import GHC.Exts (Int)--cis :: Floating a => a -> a -> Complex a-cis k n = C.cis (2 * pi * k / n)-{-# INLINE cis #-}--mkComplex :: x -> x -> Complex x-mkComplex !r !i = r :+ i-{-# INLINE mkComplex #-}+import Data.Bool (Bool,otherwise)+import Data.Bits (shiftR,shiftL,(.&.),(.|.))+import Data.Semiring (negate,(+),(*),(-))+import Control.Applicative (pure)+import Control.Monad (when)+import Data.Eq (Eq(..))+import Data.Ord (Ord(..))+import Data.Function (($),(.))+import Data.Complex (Complex(..),conjugate)+import GHC.Exts+import GHC.Real ((/))+import Control.Monad.ST (ST,runST)+import Data.Primitive.Contiguous (Contiguous,Element,Mutable)+import qualified Data.Primitive.Contiguous as Contiguous -dftMutable :: forall arr x s. (RealFloat x, Contiguous arr, Element arr (Complex x))-     => Mutable arr s (Complex x)-     -> ST s (Mutable arr s (Complex x))-dftMutable !mut = do-  !sz <- sizeMutable mut-      -  let getII !ix = (ix + sz `Prelude.div` 2) `Prelude.mod` sz -      go :: Int -- ^ i value-         -> Int -- ^ j value-         -> Complex x -- ^ accumulator-         -> ST s ()-      go !i !j !acc = if i == sz then return () else if j < sz-        then do-          let !jj = getII j-          atJJ@(r :+ _) <- read mut jj-          let real, imag, same :: x-              !same = (-2) * pi * (fromIntegral (i * j)) / (fromIntegral sz)-              !real = r * cos same-              !imag = r * sin same-              !val  = acc + mkComplex real imag-          go i (j + 1) val-        else do-          let !ii = getII i-          !_ <- write mut ii acc :: ST s ()-          go (i + 1) 0 0+{-# RULES +"fft/ifft" forall x. fft (ifft x) = x+"ifft/fft" forall x. ifft (fft x) = x+  #-} -  !_ <- go 0 0 0+-- | Radix-2 decimation-in-time fast Fourier Transform.+--   The given array must have a length that is a power of two.+fft :: forall arr. (Contiguous arr, Element arr (Complex Double))+  => arr (Complex Double)+  -> arr (Complex Double)+{-# inlinable [1] fft #-}+fft arr = if arrOK arr+  then runST $ do {+      marr <- copyWhole arr+    ; mfft marr+    ; Contiguous.unsafeFreeze marr+  }+  else Prelude.error "Data.Primitive.Contiguous.FFT.fft: bad array length" -  return mut+-- | Inverse fast Fourier transform.+ifft :: forall arr. (Contiguous arr, Element arr (Complex Double))+  => arr (Complex Double)+  -> arr (Complex Double)+{-# inlinable [1] ifft #-}+ifft arr = if arrOK arr+  then+    let lenComplex = intToComplexDouble (Contiguous.size arr)+    in cmap ((/lenComplex) . conjugate) . fft . cmap conjugate $ arr+  else Prelude.error "Data.Primitive.Contiguous.FFT.ifft: bad vector length" -dft :: forall arr x. (RealFloat x, Contiguous arr, Element arr x, Element arr (Complex x))-     => arr x-     -> arr (Complex x)-dft !a = runST $ dftInternal a+copyWhole :: forall arr s a. (Contiguous arr, Element arr a)+  => arr a+  -> ST s (Mutable arr s a)+{-# inline copyWhole #-}+copyWhole arr = do+  let len = Contiguous.size arr+  marr <- Contiguous.new len+  Contiguous.copy marr 0 arr 0 len+  pure marr --- | not in-place, also very inefficient. currently /O(n^2)/-dftInternal :: forall arr x s. (RealFloat x, Contiguous arr, Element arr x, Element arr (Complex x))-  => arr x-  -> ST s (arr (Complex x))-dftInternal !a = do-  let !sz = size a-      getII !ix = (ix + sz `Prelude.div` 2) `Prelude.mod` sz -  -  !mut <- new sz :: ST s (Mutable arr s (Complex x))+arrOK :: forall arr a. (Contiguous arr, Element arr a)+  => arr a+  -> Bool+{-# inline arrOK #-}+arrOK arr =+  let n = Contiguous.size arr+  in (1 `shiftL` log2 n) == n  -  let go :: Int -- ^ i value-         -> Int -- ^ j value-         -> Complex x -- ^ accumulator-         -> ST s ()-      go !i !j !acc = if i == sz then return () else if j < sz-        then do-          let !jj = getII j-              !atJJ = index a jj-              real, imag, same :: x-              !same = (-2) * pi * (fromIntegral (i * j)) / (fromIntegral sz)-              !real = atJJ * cos same-              !imag = atJJ * sin same-              !val  = acc + mkComplex real imag-          go i (j + 1) val-        else do-          let !ii = getII i-          !_ <- write mut ii acc :: ST s ()-          go (i + 1) 0 0--  !_ <- go 0 0 0+-- array length must be power of two. This property is not checked +mfft :: forall arr s. (Contiguous arr, Element arr (Complex Double))+  => Mutable arr s (Complex Double)+  -> ST s ()+mfft mut = do {+    len <- Contiguous.sizeMutable 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 <- Contiguous.read mut i1+                                ; let !c = Prelude.cos a+                                      !s = Prelude.sin a+                                      d = (c*xi1 - s*yi1) :+ (s*xi1 + c*yi1)+                                ; ci <- Contiguous.read mut i+                                ; Contiguous.write mut i1 (ci - d)+                                ; Contiguous.write mut i (ci + d)+                                ; butterfly (i + l2)+                              }+                      ; butterfly j+                    }+            ; flight 0 0+         }+  ; bitReverse 0 0+} -  unsafeFreeze mut+-- 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 #-} -idft :: forall arr x. (RealFloat x, Contiguous arr, Element arr x, Element arr (Complex x))-  => arr (Complex x)-  -> arr x-idft !a = runST $ idftInternal a+log2 :: Int -> Int+log2 v0 = if v0 <= 0+  then Prelude.error $ "Data.Primitive.Contiguous.FFT: nonpositive input, got " Prelude.++ Prelude.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 --- | not in-place, also very inefficient. currently /O(n^2)/-idftInternal :: forall arr x s. (RealFloat x, Contiguous arr, Element arr x, Element arr (Complex x))-  => arr (Complex x)-  -> ST s (arr x)-idftInternal !a = do-  let !sz = size a-      getII !ix = (ix + sz `Prelude.div` 2) `Prelude.mod` sz+swap :: forall arr s x. (Contiguous arr, Element arr x)+  => Mutable arr s x+  -> Int+  -> Int+  -> ST s ()+{-# inline swap #-}+swap mut i j = do+  atI <- Contiguous.read mut i+  atJ <- Contiguous.read mut j+  Contiguous.write mut i atJ+  Contiguous.write mut j atI -  !mut <- new sz :: ST s (Mutable arr s x)-  -  let go :: Int-         -> Int-         -> x-         -> ST s ()-      go !i !j !acc = if i == sz then return () else if j < sz-        then do-          let !jj = getII j-              !atJJ@(real :+ imag) = index a jj-              !sCount = fromIntegral sz-              !same = (-2) * pi * (fromIntegral (i * j)) / sCount-              !val = (real * cos same + imag * sin same) / sCount-          go i (j + 1) val-        else do-          let !ii = getII i -          !_ <- write mut ii acc :: ST s ()-          go (i + 1) 0 0+twoPi :: Double+{-# inline twoPi #-}+twoPi = 6.283185307179586 -  !_ <- go 0 0 0+intToDouble :: Int -> Double+{-# inline intToDouble #-}+intToDouble = Prelude.fromIntegral -  unsafeFreeze mut+wordToInt :: Word -> Int+{-# inline wordToInt #-}+wordToInt = Prelude.fromIntegral --- | Given a signal size, previous window, transform of previous window, and the newest value,---   compute the transform of the new window (which is just a shifted version of the previous window)---   in /O(n)/ time, in-place-overlapDFT :: forall arr x s. (RealFloat x, Contiguous arr, Element arr x, Element arr (Complex x))-  => Int   -- ^ N, signal size-  -> Mutable arr s (Complex x) -- ^ x1, original window-  -> Complex x -- ^ newest complex value-  -> Mutable arr s (Complex x) -- ^ f1, previous transform-  -> ST s (Mutable arr s (Complex x)) -- ^ f2, new transform-overlapDFT n x1 x2_N_1 f1 = do-  let !sz = fromIntegral n :: x+intToComplexDouble :: Int -> Complex Double+{-# inline intToComplexDouble #-}+intToComplexDouble = Prelude.fromIntegral -  !l <- sizeMutable f1-  !x1_0 <- read x1 0 :: ST s (Complex x)- -  let go :: Int -> ST s ()-      go !ix = if ix < l-        then do-          f1_k <- read f1 ix-          let foo' = cis (fromIntegral ix) sz-              res  = f1_k + x2_N_1 + x1_0-              fin  = foo' * res-          !_ <- write f1 ix fin-          go (ix + 1)-        else return ()-  go 0-  return f1-  +cmap :: (Contiguous arr, Element arr (Complex Double))+  => (Complex Double -> Complex Double)+  -> arr (Complex Double)+  -> arr (Complex Double)+{-# inline cmap #-}+cmap = Contiguous.map