streaming-fft (empty) → 0.1.0.0
raw patch · 7 files changed
+404/−0 lines, 7 filesdep +basedep +contiguous-fftdep +ghc-primsetup-changed
Dependencies added: base, contiguous-fft, ghc-prim, prim-instances, primitive, streaming
Files
- ChangeLog.md +5/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- src/Streaming/FFT.hs +149/−0
- src/Streaming/FFT/Internal.hs +151/−0
- src/Streaming/FFT/Types.hs +31/−0
- streaming-fft.cabal +36/−0
+ ChangeLog.md view
@@ -0,0 +1,5 @@+# Revision history for streaming-fft++## 0.1.0.0 -- YYYY-mm-dd++* First version. Released on an unsuspecting world.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2018, chessai++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of chessai nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Streaming/FFT.hs view
@@ -0,0 +1,149 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -Wall #-}++module Streaming.FFT+ ( streamFFT+ ) where++import Prelude+ ( RealFloat+ )++import Control.Monad (Monad(return))+import Control.Monad.Primitive+import Data.Complex (Complex(..))+import Data.Either (Either(..))+import Data.Eq (Eq((==)))+import Data.Function (($))+import Data.Ord (Ord(..))+import Data.Primitive.PrimArray+import Data.Primitive.Types+import GHC.Classes (modInt#)+import GHC.Num (Num(..))+import GHC.Real (fromIntegral, RealFrac(..))+import GHC.Types (Int(..))+import Streaming.FFT.Internal (initialDFT, subDFT, updateWindow', rToComplex)+import Streaming.FFT.Types (Window(..), Transform(..), Signal(..), Bin(..))+import Streaming+import Streaming.Prelude (next, yield)++data Depleted+ = NotDepleted -- ^ bin is not depleted+ | Past !Int -- ^ how many bins we have past++binDepleted :: forall e. (Num e, Ord e, RealFrac e)+ => Bin e+ -> e+ -> e+ -> Depleted+binDepleted (Bin binSize) old new =+ let !k = new - (old + fromIntegral binSize)+ in if k > 0+ then Past (floor k)+ else NotDepleted++-- [NOTE]: A drawback of the dense-stream optimisation+-- is that we must keep track of the number of bins that+-- we ingest that are 0. if too many are 0 w.r.t. the signal+-- size, then we must fall back to the /O(n log n) computation+-- until we reach another dense area of the stream. This amounts+-- to keeping an Int around that counts the number of bins that+-- were equal to zero, it gets incremented after each bin is finished+-- loading. So, there should realy be two 'thereafter' functions,+-- and 'loadInitial' should do some additional checks.+-- This is currently not the case.+loadInitial :: forall m e b. (Prim e, PrimMonad m, RealFloat e)+ => MutablePrimArray (PrimState m) (Complex e) -- ^ array to which we should allocate+ -> Bin e -- ^ bin size+ -> Signal e -- ^ signal size+ -> Int -- ^ index+ -> Int -- ^ bin accumulator+ -> e -- ^ bin pivot+ -> Int -- ^ have we finished consuming the signal+ -> Stream (Of e) m b -- first part of stream+ -> m (Stream (Of e) m b) -- stream minus original signal+loadInitial !mpa !b s@(Signal !sigSize) !ix !binAccum !binFirst !untilSig st = if (untilSig >= sigSize) then return st else do+ e <- next st+ case e of+ Left _ -> return st+ Right (x, rest) -> if ix == 0+ then loadInitial mpa b s (ix + 1) binAccum x untilSig st+ else do+ let isDepleted = binDepleted b binFirst x + case isDepleted of+ NotDepleted -> loadInitial mpa b s ix (binAccum + 1) binFirst untilSig rest+ Past i -> do+ let !k = rToComplex (fromIntegral binAccum) :: Complex e+ !_ <- writePrimArray mpa (unsafeMod (ix - 1 + untilSig) sigSize) k :: m ()+ loadInitial mpa b s (ix + i) 0 x (untilSig + 1) rest++thereafter :: forall m e b c. (Prim e, PrimMonad m, RealFloat e)+ => (Transform m e -> m c) -- ^ extract+ -> Bin e -- ^ bin size+ -> Signal e -- ^ signal size+ -> Int -- ^ index+ -> Int -- ^ have we filled a bin+ -> e -- ^ first thing in the bin+ -> Window m e -- ^ window+ -> Transform m e -- ^ transform+ -> Stream (Of e) m b+ -> Stream (Of c) m b+thereafter extract !b !s !ix !binAccum !binFirst win trans st = do+ e <- lift $ next st+ case e of+ Left r -> return r+ Right (x, rest) -> if ix == 0+ then thereafter extract b s (ix + 1) binAccum x win trans st+ else do+ let isDepleted = binDepleted b binFirst x+ case isDepleted of+ NotDepleted -> thereafter extract b s ix (binAccum + 1) binFirst win trans rest+ Past i -> do+ let k :: Complex e+ !k = rToComplex (fromIntegral binAccum)+ !trans' <- lift $ subDFT s win k trans+ !info <- lift $ extract trans'+ yield info+ -- a problem is that if too many empty bins pass,+ -- the optimised streaming-fft algorithm fails, and we+ -- need to revert (temporarily) to the original O(n log n)+ -- algorithm.+ !_ <- lift $ updateWindow' win k i+ thereafter extract b s (ix + i) 0 x win trans' rest++{-# INLINABLE streamFFT #-}+streamFFT :: forall m e b c. (Prim e, PrimMonad m, RealFloat e)+ => (Transform m e -> m c) -- ^ extraction method+ -> Bin e -- ^ bin size+ -> Signal e -- ^ signal size+ -> Stream (Of e) m b -- ^ input stream+ -> Stream (Of c) m b -- ^ output stream+streamFFT extract b s@(Signal sigSize) strm = do+ -- Allocate the one array + mpaW :: MutablePrimArray (PrimState m) (Complex e) <- lift $ newPrimArray sigSize+ let win = Window mpaW+ + -- Grab the first signal from the stream+ subStrm :: Stream (Of e) m b <- lift $ loadInitial mpaW b s 0 0 0 0 strm+ + -- Compute the transform on the signal we just grabbed+ -- so we can perform our dense-stream optimisation+ !initialT <- lift $ initialDFT win++ -- Extract information from that transform+ !initialInfo <- lift $ extract initialT+ + -- Yield that information to the new stream+ !_ <- yield initialInfo++ -- Now go+ thereafter extract b s 0 0 0 win initialT subStrm++-- | Only safe when the second argument is not 0+unsafeMod :: Int -> Int -> Int+unsafeMod (I# x#) (I# y#) = I# (modInt# x# y#)+{-# INLINE unsafeMod #-} -- this should happen anyway. trust but verify.+
+ src/Streaming/FFT/Internal.hs view
@@ -0,0 +1,151 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples #-}++{-# OPTIONS_GHC -Wall -fwarn-redundant-constraints #-}++module Streaming.FFT.Internal+ ( initialDFT+ , subDFT++ -- * some stuff (???)+ , rToComplex+ , iToComplex + , mkComplex+ , getX, getY+ + + , updateWindow + , updateWindow'+ ) where++import Control.Applicative (Applicative(pure))+import Control.Monad.Primitive+import Data.Complex hiding (cis)+import Data.Function (($))+import Data.Functor (Functor(fmap))+import Data.Primitive.PrimArray+import Data.Primitive.Types+import GHC.Num (Num(..))+import GHC.Real+import GHC.Types (Int(..))+import Prelude ()+import Data.Primitive.Instances ()+import Streaming.FFT.Types+import qualified Data.Complex as C+import qualified Data.Primitive.Contiguous.FFT as CF+import qualified Prelude as P++cis :: P.Floating e+ => e+ -> e+ -> Complex e+cis k n = C.cis (2 * P.pi * k / n)+{-# INLINE cis #-}++getX :: Complex e -> e+getX (x :+ _) = x+{-# INLINE getX #-}++getY :: Complex e -> e+getY (_ :+ y) = y+{-# INLINE getY #-}++mkComplex :: e+ -> e+ -> Complex e+mkComplex x y = x :+ y+{-# INLINE mkComplex #-}++rToComplex :: P.Num e+ => e+ -> Complex e+rToComplex e = e :+ 0+{-# INLINE rToComplex #-}++iToComplex :: P.Num e+ => e+ -> Complex e+iToComplex e = 0 :+ e+{-# INLINE iToComplex #-}++initialDFT :: forall m e. (P.RealFloat e, Prim e, PrimMonad m)+ => Window m e+ -> m (Transform m e)+initialDFT (Window !w) = fmap Transform $ stToPrim $ CF.dftMutable w+{-# INLINE initialDFT #-}++-- | Compute FFT, F2, of a Window x2 given a new sample+-- and the Transform of the old sample x1,+-- +-- IN-PLACE. (F2 is a mutated F1)+--+-- /O(n)/+subDFT :: forall m e. (P.RealFloat e, Prim e, PrimMonad m)+ => Signal e -- N+ -> Window m e -- x1+ -> Complex e+ -> Transform m e -- F1, given+ -> m (Transform m e) -- F2+subDFT (Signal n) (Window x1) x2_N_1 (Transform f1) = do+ let sz = P.fromIntegral n :: e+ l = sizeofMutablePrimArray f1 + x1_0 <- readPrimArray x1 0 :: m (Complex e)+ let go :: Int -> m ()+ go ix = if (ix P.< l)+ then do+ f1_k <- readPrimArray f1 ix+ let foo' = cis (P.fromIntegral ix) sz+ res = f1_k + x2_N_1 + x1_0+ fin = foo' * res+ writePrimArray f1 ix fin+ go (ix + 1)+ else pure ()+ go 0+ pure $ Transform f1++updateWindow' :: forall m e. (Prim e, PrimMonad m, P.RealFloat e)+ => Window m e+ -> Complex e+ -> Int -- ^ how many zeroed bins. for dense enough streams, this will be 0 most of the time+ -> m ()+updateWindow' (Window !mpa) !c !i = do+ let !sz = sizeofMutablePrimArray mpa+ !szm1 = sz - 1+ go :: Int -> m ()+ go !ix = if ix P.== szm1+ then do+ !_ <- writePrimArray mpa ix c+ P.return ()+ else if ix P.< szm1 P.&& (ix P.> szm1 P.- i)+ then do+ !_ <- writePrimArray mpa ix 0+ go (ix + 1)+ else do+ !x <- readPrimArray mpa ix+ !_ <- writePrimArray mpa (ix - 1) x+ go (ix + 1)+ go 1++updateWindow :: forall m e. (Prim e, PrimMonad m)+ => Window m e+ -> Complex e+ -> m ()+updateWindow (Window mpa1) c = do+ let !sz = sizeofMutablePrimArray mpa1+ !szm1 = sz - 1+ go :: Int -> m ()+ go !ix = if ix P.== szm1+ then do+ !_ <- writePrimArray mpa1 ix c+ P.return ()+ else do+ !x <- readPrimArray mpa1 ix+ !_ <- writePrimArray mpa1 (ix - 1) x+ go (ix + 1)+ go 1
+ src/Streaming/FFT/Types.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE GADTs #-}++{-# OPTIONS_GHC -Wall #-}++module Streaming.FFT.Types+ ( -- * types+ Signal(..)+ , Shift(..)+ , Bin(..)+ , Transform(..)+ , Window(..)+ ) where++import Control.Monad.Primitive+import Data.Complex+import Data.Primitive.PrimArray+import Prelude hiding (undefined, Rational)++-- {-# WARNING undefined "'undefined' remains in code" #-}+-- undefined :: a+-- undefined = error "Prelude.undefined"++newtype Window m e = Window+ { getWindow :: MutablePrimArray (PrimState m) (Complex e) }++newtype Transform m e = Transform+ { getTransform :: MutablePrimArray (PrimState m) (Complex e) }++newtype Signal e = Signal Int+newtype Shift e = Shift Int+newtype Bin e = Bin Int
+ streaming-fft.cabal view
@@ -0,0 +1,36 @@+name: streaming-fft+version: 0.1.0.0+synopsis: online streaming fft+description:+ online (in input and output) streaming fft algorithm+ that uses a dense-stream optimisation to reduce work+ from /O(n log n)/ to /O(n)/.+homepage: https://github.com/chessai/streaming-fft+license: BSD3+license-file: LICENSE+author: chessai+maintainer: chessai1996@gmail.com+category: Data+build-type: Simple+extra-source-files: ChangeLog.md+cabal-version: >=1.10++library+ exposed-modules:+ Streaming.FFT+ Streaming.FFT.Internal+ Streaming.FFT.Types+ build-depends:+ base >=4.9 && <5.0+ , contiguous-fft+ , ghc-prim + , prim-instances+ , primitive +-- , primitive-checked+ , streaming+ hs-source-dirs:+ src+ default-language:+ Haskell2010+ ghc-options:+ -O2