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