-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.SlidingFFT
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Sliding FFT Algorithm
--
-----------------------------------------------------------------------------

module Numeric.Transform.Fourier.SlidingFFT (sfft) where

import Data.Array
import Data.Complex

import Numeric.Transform.Fourier.FFT

-- Sliding FFT algorithm.  We assume that the head of the list is the
-- oldest sample, and the last element is the newest sample.  This is why
-- we need the reverse.  By doing this we can abstract things like A/D
-- converters as infinite lists.

-- The only published reference I have seen for this is the TI TMS320C3x
-- General-Purpose Applications (SPRU194).  You can also check out
-- comp.dsp.  The author, Keith Larson, hangs out there.

-- The type of (!!) forces the type signatures to use Int instead of
-- (Integral a)

{-# specialize sfft :: Int -> [Complex Float] -> [Array Int (Complex Float)] #-}
{-# specialize sfft :: Int -> [Complex Double] -> [Array Int (Complex Double)] #-}

-- | Sliding FFT

sfft :: RealFloat a => Int -- ^ N
     -> [Complex a] -- ^ x[n]
     -> [Array Int (Complex a)] -- ^ [X[k]]

sfft :: forall a.
RealFloat a =>
Int -> [Complex a] -> [Array Int (Complex a)]
sfft Int
_ [] = forall a. HasCallStack => [Char] -> a
error [Char]
"sfft: input must have at least on value"
sfft Int
n (Complex a
x:[Complex a]
xs) = Array Int (Complex a)
x' forall a. a -> [a] -> [a]
: forall a.
RealFloat a =>
Int
-> Complex a
-> [Complex a]
-> Array Int (Complex a)
-> [Array Int (Complex a)]
sfft' Int
n Complex a
x [Complex a]
xs Array Int (Complex a)
x'
    where x' :: Array Int (Complex a)
x' = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b)
fft forall a b. (a -> b) -> a -> b
$ forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
nforall a. Num a => a -> a -> a
-Int
1) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take Int
n (Complex a
xforall a. a -> [a] -> [a]
:[Complex a]
xs)

{-# specialize sfft' :: Int -> Complex Float -> [Complex Float] -> Array Int (Complex Float) -> [Array Int (Complex Float)] #-}
{-# specialize sfft' :: Int -> Complex Double -> [Complex Double] -> Array Int (Complex Double) -> [Array Int (Complex Double)] #-}

sfft' :: RealFloat a => Int
     -> Complex a
     -> [Complex a]
     -> Array Int (Complex a)
     -> [Array Int (Complex a)]
sfft' :: forall a.
RealFloat a =>
Int
-> Complex a
-> [Complex a]
-> Array Int (Complex a)
-> [Array Int (Complex a)]
sfft' Int
n Complex a
xn (Complex a
x:[Complex a]
xs)  Array Int (Complex a)
x' | forall a. Int -> [a] -> Bool
enough Int
n (Complex a
xforall a. a -> [a] -> [a]
:[Complex a]
xs) = Array Int (Complex a)
x'' forall a. a -> [a] -> [a]
: forall a.
RealFloat a =>
Int
-> Complex a
-> [Complex a]
-> Array Int (Complex a)
-> [Array Int (Complex a)]
sfft' Int
n Complex a
x [Complex a]
xs Array Int (Complex a)
x''
		      | Bool
otherwise       = []
    where x'' :: Array Int (Complex a)
x'' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
nforall a. Num a => a -> a -> a
-Int
1) [ Complex a
x0 forall a. Num a => a -> a -> a
- Complex a
xn forall a. Num a => a -> a -> a
+ Array Int (Complex a)
x'forall i e. Ix i => Array i e -> i -> e
!Int
i forall a. Num a => a -> a -> a
* forall {a} {a}. (Floating a, Integral a) => a -> Complex a
w Int
i | Int
i <- [Int
0..(Int
nforall a. Num a => a -> a -> a
-Int
1)] ]
          x0 :: Complex a
x0  = [Complex a]
xs forall a. [a] -> Int -> a
!! (Int
nforall a. Num a => a -> a -> a
-Int
2)
	  w :: a -> Complex a
w a
i = forall a. Floating a => a -> Complex a
cis forall a b. (a -> b) -> a -> b
$ -a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
sfft' Int
_ Complex a
_ [] Array Int (Complex a)
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"sfft': input must have at least on value"

-- We can't use Prelude.length because we may be operating on infinite,
-- or ginormous lists.  So enough will return True is there is enough
-- data to perform the next FFT update, or False if there is not enough.

enough :: Int -> [a] -> Bool
enough :: forall a. Int -> [a] -> Bool
enough Int
n [a]
xs  =  Int
nforall a. Ord a => a -> a -> Bool
<=Int
0 Bool -> Bool -> Bool
|| Bool -> Bool
not (forall (t :: * -> *) a. Foldable t => t a -> Bool
null (forall a. Int -> [a] -> [a]
drop (Int
nforall a. Num a => a -> a -> a
-Int
1) [a]
xs))

{-
alternative:

enough 0 _      = True
enough _ []     = False
enough n (_:xs) = enough (n-1) xs
-}