HaskellでFFT
HaskellでstockhamのFFTを作ってみた。
Listなので、64Kサンプルになると遅くなる。ここから、どこまで速くなるか試みてみる。
module Numerical.FFT ( fft1d ) where import Debug.Trace import System.Environment import System.Console.GetOpt import System.IO import Data.Maybe import Data.Complex import Data.Bits {- 時間間引きのstockham FFT -} {- スキップ数に応じたListを作ります -} sliceTwoLineList :: Int -> [[Complex a]] -> [[[Complex a]]] sliceTwoLineList n tl = sliceTwoLineListIter tl [] where sliceTwoLineListIter :: [[Complex a]]-> [[[Complex a]]] -> [[[Complex a]]] sliceTwoLineListIter tl ltl | null (tl!!0) == True = ltl | otherwise = sliceTwoLineListIter [drop n (tl!!0),drop n (tl!!1)] (ltl ++ [[take n (tl!!0),take n (tl!!1)]]) {- FFTの本体 -} stockham :: (Show a,RealFloat a) => [Complex a] -> [Complex a] stockham inDataList = stockhamIter (div (length inDataList) 2) 1 inDataList where stockhamIter :: (Show a,RealFloat a) => Int -> Int -> [Complex a] -> [Complex a] stockhamIter n skip inDataList | n == 1 = butterflyCalc skip butterflySrc | otherwise = stockhamIter (div n 2) (skip*2) (butterflyCalc skip butterflySrc) where splitedDataList = let s = splitAt (div (length inDataList) 2) inDataList in [fst s,snd s] {- データを前半後半に分けて、スキップ数に分けます -} butterflySrc = sliceTwoLineList skip splitedDataList butterflyCalc :: (Show a,RealFloat a) => Int -> [[[Complex a]]] -> [Complex a] butterflyCalc skip src = concat $ map (\x -> bufferflyCalc' x) src where multipleTwiddle :: (Show a, RealFloat a) => (Complex a -> Complex a -> Complex a) -> Int -> Complex a -> Complex a -> Complex a multipleTwiddle f m op1 op2 = let w = cis ((-pi*(fromIntegral m))/(fromIntegral skip)) -- e^(-(2*pi*m)/(2*skip)) in f op1 (op2 * w) bufferflyCalc' :: (Show a, RealFloat a) => [[Complex a]] -> [Complex a] {- 前半は加算、後半は減算 -} bufferflyCalc' src = (bufferflyCalcIter (multipleTwiddle (+)) 0 src []) ++ (bufferflyCalcIter (multipleTwiddle (-)) 0 src []) where bufferflyCalcIter :: (Show a, RealFloat a) => (Int -> Complex a -> Complex a -> Complex a) -> Int -> [[Complex a]] -> [Complex a] -> [Complex a] bufferflyCalcIter f n src l | null (src!!0) == True = l | otherwise = bufferflyCalcIter f (n+1) [drop 1 (src!!0),drop 1 (src!!1)] (l ++ [f n (head (src!!0)) (head (src!!1))]) {- 波形データを引き取って周波数成分を返す -} fft1dCore :: (Show a, RealFloat a) => [Complex a] -> Int -> Bool -> [Complex a] fft1dCore inData numOfTotal bInverse = if bInverse then let inData' = map conjugate inData convData = stockham inData in map (/ ((fromIntegral numOfTotal) :+ 0) ) $ map conjugate convData else stockham inData fft1d :: (Show a, RealFloat a) => [Complex a] -> Bool -> Maybe [Complex a] fft1d inData bInverse = let n = length inData in if (n .&. (n-1) == 0) then Just (fft1dCore inData n bInverse) else Nothing -- Dataの長さが2のべき乗でありません