diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,9 @@ # Revision history for streaming-fft -## 0.1.0.0 -- YYYY-mm-dd +## 0.1.0.1 -- 2018-10-24 + +* Documentation improvements + +## 0.1.0.0 -- 2018-10-07 * First version. Released on an unsuspecting world. diff --git a/src/Streaming/FFT.hs b/src/Streaming/FFT.hs --- a/src/Streaming/FFT.hs +++ b/src/Streaming/FFT.hs @@ -5,7 +5,12 @@ {-# OPTIONS_GHC -Wall #-} module Streaming.FFT - ( streamFFT + ( -- * streaming fft + streamFFT + -- * types + , Transform(..) + , Bin(..) + , Signal(..) ) where import Prelude @@ -114,20 +119,79 @@ !_ <- 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 +-- | 'streamFFT' is based off ideas from signal processing, with an optimisation +-- outlined in . +-- Here, I will give you an outline of how this works. The idea is that we +-- have a stream of data, which we will divide into 'Signal's, and each 'Signal' +-- is something for which we want to compute the DFT. Each signal is divided into +-- 'Bin's (more on this later, but you can just think of 'Bin's as a chunk of a +-- 'Signal', where all the chunks are of equal length). We treat our stream not as +-- contiguous blocks of 'Signal's, but as overlapping 'Signal's, where each overlap +-- is one 'Bin'-length. The motivation for the blog post is to reduce the work of +-- this overlap; they show a way to compute the DFT of each 'Signal' subsequent +-- to the initial in /O(n)/ time, instead of the typical /O(n log n)/ time, +-- by abusing the overlap. +-- +-- Consider you would like to compute the Fourier Transform of the signal +-- +-- \[ +-- x_{i-n+1}, x_{i-n+2}, ..., x_{i-1}, x_{i}. +-- \] +-- +-- However this means that when you receive \( x_{i+1} \), you'll be the computing +-- the Fourier Transform of +-- +-- \[ +-- x_{i-n+2}, x_{i-n+3}, ..., x_{i}, x_{i+1}, +-- \] +-- +-- which is almost identical to the first sequence. How do we avoid extra work? +-- +-- Assume data windows to be of length \( N \) (this corresponds to the number of +-- 'Bin's in the 'Signal'). Let +-- the original data window be \( x_{1} \), whose first sample is \( x_{old} = x_{1}[0] \). +-- (here, \( a[k] \) is used to denote accessing the \( (k-1)th \) element from +-- a sequence \( a \) ). Let your new data window be denoted as \( x_{2} \), whose +-- bins are one left-shifted version of \( x_{1} \), i.e. +-- \( x_{2}[k] = x_{1}[k+1]\) for \(k = 0, 1, ... N - 2 \), plus a new arrived datum to +-- position \( N - 1 \), which is denoted as \( x_{new} = x_{2}[N - 1]\). +-- +-- The following will compute the N-point DFT, \( X_{2} \) of the new data set +-- \( x_{2} \) from that of the already computed and stored N-point DFT +-- \( X_{1} \) of the old data set \( x_{1} \): +-- +-- \[ +-- X{2}[k] = e^{2 \pi i k / N} * (X{1}[k] + (x_{new} - x_{old})) +-- \] +-- +-- for each \( k = 0, 1, ..., N - 1 \). This updated computation of \( X{2} \) +-- pre-computed \( X{1} \) requires \( N \) complex multiplications and \( N \) +-- real additions. Compared to a direct N-point DFT which requires \( N log_{2}(N) \) +-- complex multiply-accumulate operations, this is an improvement by a factor of +-- \( log_{2}(N) \), which for example at N=1024 would translate to a speedup of +-- about 10. +-- +-- Another advantage of this algorithm as this it is amenable to being done in-place. +-- `streamFFT` in fact does do this, and for that reason allocations are kept to an +-- absolute minimum. +-- +-- +streamFFT :: forall m a b c. (Prim a, PrimMonad m, RealFloat a) + => (Transform m a -> m c) -- ^ extraction method. This is a function that takes a 'Transform' + -- and produces (or 'extracts') some value from it. It is used + -- to produce the values in the output stream. + -> Bin a -- ^ bin size + -> Signal a -- ^ signal size + -> Stream (Of a) m b -- ^ input stream -> Stream (Of c) m b -- ^ output stream +{-# INLINABLE streamFFT #-} streamFFT extract b s@(Signal sigSize) strm = do -- Allocate the one array - mpaW :: MutablePrimArray (PrimState m) (Complex e) <- lift $ newPrimArray sigSize + mpaW <- 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 + subStrm <- 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 @@ -146,4 +210,3 @@ unsafeMod :: Int -> Int -> Int unsafeMod (I# x#) (I# y#) = I# (modInt# x# y#) {-# INLINE unsafeMod #-} -- this should happen anyway. trust but verify. - diff --git a/src/Streaming/FFT/Types.hs b/src/Streaming/FFT/Types.hs --- a/src/Streaming/FFT/Types.hs +++ b/src/Streaming/FFT/Types.hs @@ -1,31 +1,33 @@ -{-# LANGUAGE GADTs #-} +{-# LANGUAGE GADTSyntax #-} +{-# LANGUAGE TypeInType #-} {-# OPTIONS_GHC -Wall #-} module Streaming.FFT.Types ( -- * types Signal(..) - , Shift(..) , Bin(..) , Transform(..) , Window(..) ) where +import Data.Kind (Type) 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) } +-- | A 'Window' is a mutable primitive array of 'Complex' values, +-- over which we compute the DFT. +newtype Window :: (Type -> Type) -> Type -> Type where + Window :: MutablePrimArray (PrimState m) (Complex e) -> Window m e -newtype Transform m e = Transform - { getTransform :: MutablePrimArray (PrimState m) (Complex e) } +-- | A 'Transform' is a Mutable primitive array of 'Complex' values, +-- the result of taking the DFT of a 'Window'. +newtype Transform :: (Type -> Type) -> Type -> Type where + Transform :: MutablePrimArray (PrimState m) (Complex e) -> Transform m e +-- | Your signal size. newtype Signal e = Signal Int -newtype Shift e = Shift Int +-- | Your bin size. newtype Bin e = Bin Int diff --git a/streaming-fft.cabal b/streaming-fft.cabal --- a/streaming-fft.cabal +++ b/streaming-fft.cabal @@ -1,5 +1,5 @@ name: streaming-fft -version: 0.1.0.0 +version: 0.1.0.1 synopsis: online streaming fft description: online (in input and output) streaming fft algorithm