HaskellでFFT.Vector編
二次元配列を一次元で表現するために、二重ループが入るが、その際のVectorのappendで時間がかかる…
module Numerical.FFT ( fft1d ) where import Prelude hiding ((++),length,map) import Data.Vector import Data.Maybe import Data.Complex import Data.Bits {- 時間間引きのstockham FFT -} doubleLoop :: (Num a) => (Int -> Int -> a) -> Int -> Int -> Vector a doubleLoop f jMax iMax = jIter 0 empty where jIter j holdData | j < jMax = jIter (j+1) (holdData ++ (iIter iMax)) | otherwise = holdData where iIter i = generate i (\i -> f i j) butterflyCalc :: (RealFloat a) => Vector (Complex a) -> Vector (Complex a) -> Vector (Complex a) butterflyCalc twiddleTable srcDatas = butterflyCalcIter (length srcDatas) 1 srcDatas where butterflyCalcIter n m srcDatas | n == 1 = srcDatas | otherwise = let n' = (div n 2) m' = m * 2 in butterflyCalcIter n' m' $ doubleLoop calc n' m' where calc i j = let n' = div n 2 l = n' * m -- = div (length srcData) 2 in if(i < m) then let w = twiddleTable!(i*n') in srcDatas!(j*m+i)+(srcDatas!(j*m+l+i))*w else let w = twiddleTable!((i-m)*n') in srcDatas!(j*m+i-m)-(srcDatas!(j*m+l+i-m))*w fft1dCore :: (RealFloat a) => Vector (Complex a) -> Bool -> Vector (Complex a) fft1dCore inDatas bInverse = let numOfSample = length inDatas signOfRootOfUnity = if bInverse then 1 else -1 twiddlesTable = generate numOfSample (\i -> cis (signOfRootOfUnity*2*pi*(fromIntegral i)/(fromIntegral numOfSample))) --回転子のωNのテーブルを作ります。 outDatas = butterflyCalc twiddlesTable inDatas in if bInverse then map (\z -> z / (fromIntegral numOfSample)) outDatas else outDatas fft1d :: (Show a, RealFloat a) => Vector (Complex a) -> Bool -> Maybe (Vector (Complex a)) fft1d inDatas bInverse = let n = length inDatas in if (n .&. (n-1) == 0) then Just (fft1dCore inDatas bInverse) else Nothing -- Dataの長さが2のべき乗でありません