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のべき乗でありません