FFT.hs

-- | This module implements the Fast Fourier Transformation purely in Haskell.
module Sound.Hommage.FFT
 ( 
 -- * FFT
   fft
 , fft'
 , fftt 

 -- * FFT
 , fftc
 , fftc' 

 , ffttv 

 , fftco 
 , fftco'

 , rect_transform 
 , rect_filter 
 , rect_filter' 

 , i 
 , w 
 , ws 
 , map2 
 , zip2 
 , evens 
 , odds 
 , drop_odds 
 , appendpair
-- , fitlastlength 
 , reorder 
 , reorder_init 
 , reorder_init'
-- , reorder_overlap 
 , butterfly 
 , butterfly'
 , fft_level 
 , fft_merge 
 , fft_init_overlap 
 , fft_level' 
 , fft_merge' 
 , fft_init 
 , fft_overlap_loop 
 , expandComplex 
 , overlap_curve 
 , mix_overlap 

 )
 where
import Data.Complex
import Sound.Hommage.Misc
-------------------------------------------------------------------------------
-- | The complex value 'i' = @(0 :+ 1)@
i :: Complex Double
i = 0.0 :+ 1.0

-- | The n-th root of 1
w :: Int -> Complex Double
w n = exp (-2 * i * pi / fromIntegral n)

-- | the 2n-th root with exponents 0, 1, .. n. False=inverse (exponents are negated)
ws :: Bool -> Int -> [Complex Double]
ws b n = take n $ map (w1^^) ix 
 where
  w1 = exp (-i * pi / fromIntegral n)
  ix | b         = [0, 1 ..]
     | otherwise = [0, (-1) ..]
-------------------------------------------------------------------------------
map2 :: (a -> a -> b) -> [a] -> [b]
map2 f (x:y:r) = f x y : map2 f r
map2 _ _       = []

map2opt :: (a -> a -> b) -> a -> [a] -> [b]
map2opt f a = loop 
 where
  loop (x:y:r) = f x y : loop r
  loop [x]     = [f x a]
  loop _       = []

zip2 :: [a] -> [a] -> [a]
zip2 (x:xs) (y:ys) = x : y : zip2 xs ys
zip2 _ _ = []

evens (x:xs) = x : odds xs
evens []     = []

odds (_:xs) = evens xs
odds []     = []

-- | returns [] if argument has zero or one element.
drop_odds (x:y:r) = x : evens r
drop_odds _       = []

appendpair :: ([a], [a]) -> [a]
appendpair (x,y) = x++y

--fitlastlength :: Int -> a -> [[a]] -> [[a]]
--fitlastlength n a = loop
-- where
--  loop (x:xs) = case loop xs of
--                [] -> [take n (x ++ repeat a)]
--                xs -> x : xs
--  loop []     = []
-------------------------------------------------------------------------------
-- | list is grouped into sublists with length N (must be power of 2) and bitwise reverse order
reorder :: Int -> [[a]] -> [[a]]
reorder 0 = id
reorder n = map2 zip2 . reorder (n-1)

reorder_r2 :: Int -> a -> [[a]] -> [[a]]
reorder_r2 n a = loop n
 where
  loop 0 = id
  loop n = map2opt zip2 (repeat a) . loop (n-1)

reorder_init (x1:x2:x3:x4:xs) = [x1 :+ x2, x3 :+ x4] : reorder_init xs
reorder_init [x1,x2,x3]       = [[x1 :+ x2, x3 :+ 0]]
reorder_init [x1,x2]          = [[x1 :+ x2, 0]]
reorder_init [x1]             = [[x1 :+ 0, 0]]
reorder_init _                = []

reorder_init' :: [a] -> [[a]]
reorder_init' (x:y:r) = [x,y] : reorder_init' r
reorder_init' _       = []

reorder_init_r2' :: a -> [a] -> [[a]]
reorder_init_r2' a = loop
 where
  loop (x:y:r) = [x,y] : loop r
  loop [x]     = [[x,a]]
  loop _       = []
  
--reorder_overlap :: [[a]] -> [a]
--reorder_overlap (xs:ys:r) = (zip2 xs ys) ++ reorder_overlap (ys:r)
--reorder_overlap _       = []
-------------------------------------------------------------------------------
butterfly :: Complex Double ->  Complex Double -> Complex Double -> (Complex Double, Complex Double)
butterfly (wr :+ wi) (xr :+ xi) (yr :+ yi) = (a, b)
 where
  wyr = wr*yr-wi*yi 
  wyi = wr*yi+yr*wi
  a = ( (xr + wyr) * 0.5 :+ (xi + wyi) * 0.5 )
  b = ( (xr - wyr) * 0.5 :+ (xi - wyi) * 0.5 )

butterfly' :: Complex Double ->  Complex Double -> Complex Double -> (Complex Double, Complex Double)
butterfly' (wr :+ wi) (xr :+ xi) (yr :+ yi) = (a, b)
 where
  wyr = wr*yr-wi*yi 
  wyi = wr*yi+yr*wi
  a = ( (xr + wyr) :+ (xi + wyi) )
  b = ( (xr - wyr) :+ (xi - wyi) )
-------------------------------------------------------------------------------
-- | FFT transformation. Input is grouped into overlapping parts of 2^(N+2) reals and mapped to sublists 
--   with 2^(N+1) complex numbers. 
fft :: Int -> [Double] -> [[Complex Double]]
fft n = fft_level (ws True (2^n)) . fft_init_overlap . reorder_r2 (n-1) 0 . reorder_init

fft_level :: [Complex Double] -> [[Complex Double]] -> [[Complex Double]] 
fft_level [_] = id
fft_level ws  = map2 (\x -> appendpair . fft_merge ws x) . fft_level (drop_odds ws)

fft_merge :: [Complex Double] -> [Complex Double] -> [Complex Double] -> ([Complex Double], [Complex Double])
fft_merge (w:ws) (x:xs) (y:ys) = let (xs', ys') = fft_merge ws xs ys
                                     (x', y') = butterfly w x y
                                 in (x' : xs', y' : ys')
fft_merge _ _ _ = ([], [])

fft_init_overlap (x:y:r) = (zipWith (\x y -> [(x + y)*0.5, (x - y)*0.5]) x y) ++ fft_init_overlap (y:r)
fft_init_overlap _       = []

fft_init (x:y:r) = (zipWith (\x y -> [(x + y)*0.5, (x - y)*0.5]) x y) ++ fft_init r
fft_init _       = []
-------------------------------------------------------------------------------
-- | Inverse FFT transformation. Complex input is grouped into parts with length 2^(N+1) and mapped to 
--   sublists with 2^(N+2) reals, which are overlapped and mixed.
fft' :: Int -> [Complex Double] -> [Double]
fft' n = mix_overlap (2^(n+2)) . map expandComplex . fft_level' (ws False (2^n)) . fft_init' . reorder (n-1) . reorder_init'

fft_level' :: [Complex Double] -> [[Complex Double]] -> [[Complex Double]] 
fft_level' [_] = id
fft_level' ws  = map2 (\x -> appendpair . fft_merge' ws x) . fft_level' (drop_odds ws)

fft_merge' :: [Complex Double] -> [Complex Double] -> [Complex Double] -> ([Complex Double], [Complex Double])
fft_merge' (w:ws) (x:xs) (y:ys) = let (xs', ys') = fft_merge' ws xs ys
                                      (x', y') = butterfly' w x y
                                  in (x' : xs', y' : ys')
fft_merge' _ _ _ = ([], [])

fft_init' (x:y:r) = (zipWith (\x y -> [x + y, x - y]) x y) ++ fft_init' r
fft_init' _       = []

fft_overlap_loop (x:r@(y:_)) = (concat $ zipWith (\x y -> [x + y, x - y]) x y) : fft_overlap_loop r
fft_overlap_loop _           = []

expandComplex ((x:+y):r) = x : y : expandComplex r
expandComplex _ = []
-------------------------------------------------------------------------------
-- | Creates a fade with length N
overlap_curve :: Int -> [Double]
overlap_curve n = take n $ map (\k -> 0.5 - 0.5 * cos (fromIntegral k * 2 * pi / fromIntegral n )) [0..]

-- | Overlaps a sequence of parts of length N (overlaps by N\/2).
mix_overlap :: Int -> [[Double]] -> [Double]
mix_overlap n xs = merge (+) l r
  where 
  c = cycle $ overlap_curve n
  l = zipWith (*) c $ concat $ evens xs
  r = replicate (div n 2) 0.0 ++ (zipWith (*) c $ concat $ odds xs)
-------------------------------------------------------------------------------
-- | FFT transformation with overlapping and windowing. 
--   Argument function maps the coefficients. Windowsize: 2^(N+2), i. e. 2^(N+1) coefficients
fftt :: Int -> ([[Complex Double]] -> [Complex Double]) -> [Double] -> [Double]
fftt n f = 
 mix_overlap (2^(n+2)) . 
 map expandComplex . 
 fft_level' (ws False (2^n)) . 
 fft_init' . 
 reorder (n-1) . 
 reorder_init' .
 f .
 fft_level (ws True (2^n)) . 
 fft_init_overlap . 
 reorder_r2 (n-1) 0 . 
 reorder_init

ffttv :: Int -> [Double] -> [Double] -> [Double]
ffttv n vs = fftt n (zipWith (\v (x:+y) -> (v*x) :+ (v*y)) vs . concat . map (\((a:+_):r) -> (a:+0):r))
-------------------------------------------------------------------------------
-- | FFT transformation for complex input (segments of length 2^n). 
--   No overlapping or windowing.
fftc :: Int -> [Complex Double] -> [[Complex Double]]
fftc n = fft_level (ws True (2^n)) . fft_init . reorder_r2 (n-1) 0 . reorder_init_r2' 0

-- | Inverse FFT transformation for complex input (segments of length 2^n). 
--   No overlapping or windowing.
fftc' :: Int -> [Complex Double] -> [[Complex Double]]
fftc' n = fft_level' (ws False (2^n)) . fft_init' . reorder (n-1) . reorder_init'
-------------------------------------------------------------------------------
-- | FFT for complex input with overlapping. Segment-size: 2^(n+1)
fftco :: Int -> [Complex Double] -> [[Complex Double]]
fftco n = fft_level (ws True (2^n)) . fft_init_overlap . reorder_r2 (n-1) 0 . reorder_init_r2' 0

fftco' :: Int -> [Complex Double] -> [Complex Double]
fftco' n = mix_overlapc (2^(n+1)) . fft_level' (ws False (2^n)) . fft_init' . reorder (n-1) . reorder_init'

-- | Overlaps a sequence of parts of length N (overlaps by N\/2).
mix_overlapc :: Int -> [[Complex Double]] -> [Complex Double]
mix_overlapc n xs = merge (+) l r
  where 
  c = cycle $ map (:+0) $ overlap_curve n
  l = zipWith (*) c $ concat $ evens xs
  r = replicate (div n 2) 0.0 ++ (zipWith (*) c $ concat $ odds xs)
-------------------------------------------------------------------------------

-- | A self-inverse transformation similar to FFT but with a simple butterfly
--   operation that uses always W=1. Modifying the data between application
--   and inverse is similar to filtering but the result will be built up from
--   rectangle waves instead of sinus waves.
rect_transform :: Floating a => Int -> [a] -> [[a]]
rect_transform n = run n . map (:[])
 where
  sq = 1 / sqrt 2
  add x y = (x+y) * sq
  sub x y = (x-y) * sq
  run 0 = id
  run n = loop . run (n-1)
  loop (x:y:r) = iter True x y : loop r
  loop [x]     = [iter True x (repeat 0)]
  loop _       = []
  iter True  (x:xs) (y:ys) = add x y : sub x y : iter False xs ys
  iter False (x:xs) (y:ys) = sub x y : add x y : iter True  xs ys
  iter _     _      _      = []

rect_filter :: Floating a => Int -> ([[a]] -> [[a]]) -> [a] -> [a]
rect_filter n f =  concat . rect_transform n . concat . f . rect_transform n

rect_filter' :: Floating a => Int -> [[a]] -> [a] -> [a]
rect_filter' n vs =  concat . rect_transform n . concat . zipWith f vs . rect_transform n
 where
  f v c = zipWith (*) (v ++ repeat 0) c