summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorchessai <>2018-10-24 15:48:00 (GMT)
committerhdiff <hdiff@hdiff.luite.com>2018-10-24 15:48:00 (GMT)
commita8c84be663b402355ccf649d32611185ce5ed828 (patch)
tree9457b90d518558411f7ff3038796dec9d83ffbaa
parent1b152f862402e29ed90bb4e07674dec7d58c8c87 (diff)
version 0.1.0.1HEAD0.1.0.1master
-rwxr-xr-xChangeLog.md6
-rw-r--r--src/Streaming/FFT.hs83
-rw-r--r--src/Streaming/FFT/Types.hs24
-rw-r--r--streaming-fft.cabal2
4 files changed, 92 insertions, 23 deletions
diff --git a/ChangeLog.md b/ChangeLog.md
index bb6bd77..ac11d0b 100755
--- 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
index 3ccc93e..a4f2cea 100644
--- 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 @@ thereafter extract !b !s !ix !binAccum !binFirst win trans st = do
!_ <- 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 <https://www.dsprelated.com/showarticle/776.php this blog post>.
+-- 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 @@ streamFFT extract b s@(Signal sigSize) strm = do
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
index f481d71..38db442 100644
--- 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"
+-- | 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 Window m e = Window
- { getWindow :: MutablePrimArray (PrimState m) (Complex 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
index 6d80296..48f78ed 100644
--- 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