HaskellでFFT.Array編

Data.Arrayを使って、HaskellでFFT.
64KのサンプルでList版から101.1sec->0.73secに短縮。
大量のサンプル数を突っ込むとギブアップする。

module Numerical.FFT
(
 fft1d
) where

import qualified Data.Vector as V 
import Data.Array
import Data.Array.IArray (amap)
import Data.Maybe
import Data.Complex
import Data.Bits

{- 時間間引きのstockham FFT -}

butterflyCalc :: (RealFloat a, Show a) => V.Vector (Complex a) -> Array (Int,Int) (Complex a) -> Array (Int,Int) (Complex a)
butterflyCalc twiddlesTable srcArray = 
  let (n',m') = snd $ bounds srcArray
      n = div (n'+1) 2
      m = (m'+1) * 2
      l = m'+1
  in array ((0,0),(n-1,m-1)) [((i,j), if (j < l)
                                      then
                                        let w = ((V.!) twiddlesTable (n*j))
                                        in srcArray!(i,j)+(srcArray!(i+n,j))*w
                                      else 
                                        let w = ((V.!) twiddlesTable (n*(j-l)))
                                        in srcArray!(i,j-l)-(srcArray!(i+n,j-l))*w) |i<-[0..n-1],j<-[0..m-1]]
 
butterflyCalcIter :: (RealFloat a,Num a,Show a) => V.Vector (Complex a) -> Array (Int,Int) (Complex a) -> Array (Int,Int) (Complex a)
butterflyCalcIter twiddlesTable currentArray 
  | ((fst $ snd $ bounds currentArray) == 0) = currentArray
  | otherwise = butterflyCalcIter twiddlesTable $ butterflyCalc twiddlesTable currentArray
    
fft1dCore :: (Show a, RealFloat a) => Array (Int,Int) (Complex a) -> Bool -> Array (Int,Int) (Complex a)
fft1dCore inDatas bInverse = 
  let numOfSample = rangeSize $ bounds inDatas
      signOfRootOfUnity = if bInverse then 1 else -1 
      twiddlesTable = V.generate numOfSample (\i -> cis (signOfRootOfUnity*2*pi*(fromIntegral i)/(fromIntegral numOfSample))) --回転子のωNのテーブルを作ります。
      outDatas = butterflyCalcIter twiddlesTable inDatas
  in if bInverse then amap (\z -> z/(fromIntegral numOfSample)) outDatas else outDatas 

fft1d :: (Show a, RealFloat a) => Array (Int,Int) (Complex a) -> Bool -> Maybe (Array (Int,Int) (Complex a))
fft1d inDatas bInverse = 
  let n = rangeSize $ bounds inDatas
  in if (n .&. (n-1) == 0)
     then Just (fft1dCore inDatas bInverse)
     else Nothing -- Dataの長さが2のべき乗でありません