module Sound.Hommage.FFT
(
fft
, fft'
, fftt
, fftc
, fftc'
, ffttv
, fftco
, fftco'
, rect_transform
, rect_filter
, rect_filter'
, i
, w
, ws
, map2
, zip2
, evens
, odds
, drop_odds
, appendpair
, reorder
, reorder_init
, reorder_init'
, 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
i :: Complex Double
i = 0.0 :+ 1.0
w :: Int -> Complex Double
w n = exp (2 * i * pi / fromIntegral n)
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 [] = []
drop_odds (x:y:r) = x : evens r
drop_odds _ = []
appendpair :: ([a], [a]) -> [a]
appendpair (x,y) = x++y
reorder :: Int -> [[a]] -> [[a]]
reorder 0 = id
reorder n = map2 zip2 . reorder (n1)
reorder_r2 :: Int -> a -> [[a]] -> [[a]]
reorder_r2 n a = loop n
where
loop 0 = id
loop n = map2opt zip2 (repeat a) . loop (n1)
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 _ = []
butterfly :: Complex Double -> Complex Double -> Complex Double -> (Complex Double, Complex Double)
butterfly (wr :+ wi) (xr :+ xi) (yr :+ yi) = (a, b)
where
wyr = wr*yrwi*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*yrwi*yi
wyi = wr*yi+yr*wi
a = ( (xr + wyr) :+ (xi + wyi) )
b = ( (xr wyr) :+ (xi wyi) )
fft :: Int -> [Double] -> [[Complex Double]]
fft n = fft_level (ws True (2^n)) . fft_init_overlap . reorder_r2 (n1) 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 _ = []
fft' :: Int -> [Complex Double] -> [Double]
fft' n = mix_overlap (2^(n+2)) . map expandComplex . fft_level' (ws False (2^n)) . fft_init' . reorder (n1) . 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 _ = []
overlap_curve :: Int -> [Double]
overlap_curve n = take n $ map (\k -> 0.5 0.5 * cos (fromIntegral k * 2 * pi / fromIntegral n )) [0..]
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)
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 (n1) .
reorder_init' .
f .
fft_level (ws True (2^n)) .
fft_init_overlap .
reorder_r2 (n1) 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))
fftc :: Int -> [Complex Double] -> [[Complex Double]]
fftc n = fft_level (ws True (2^n)) . fft_init . reorder_r2 (n1) 0 . reorder_init_r2' 0
fftc' :: Int -> [Complex Double] -> [[Complex Double]]
fftc' n = fft_level' (ws False (2^n)) . fft_init' . reorder (n1) . reorder_init'
fftco :: Int -> [Complex Double] -> [[Complex Double]]
fftco n = fft_level (ws True (2^n)) . fft_init_overlap . reorder_r2 (n1) 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 (n1) . reorder_init'
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)
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 = (xy) * sq
run 0 = id
run n = loop . run (n1)
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