----------------------------------------------------------------------------- -- | -- Module : Numeric.Transform.Fourier.FFTHard -- Copyright : (c) Matthew Donadio 2003 -- License : GPL -- -- Maintainer : m.p.donadio@ieee.org -- Stability : experimental -- Portability : portable -- -- Hard-coded FFT transforms -- ----------------------------------------------------------------------------- module Numeric.Transform.Fourier.FFTHard where import Data.Array import Data.Complex -- These are the hard coded DFT's borrowed from FFTW {-# specialize fft'2 :: Array Int (Complex Float) -> Array Int (Complex Float) #-} {-# specialize fft'2 :: Array Int (Complex Double) -> Array Int (Complex Double) #-} -- | Length 2 FFT fft'2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n] -> Array a (Complex b) -- ^ X[k] fft'2 a = array (0,1) [ (0, ((tmp1 + tmp2) :+ (tmp3 + tmp4))), (1, ((tmp1 - tmp2) :+ (tmp3 - tmp4) )) ] where tmp1 = realPart (a!0) tmp3 = imagPart (a!0) tmp2 = realPart (a!1) tmp4 = imagPart (a!1) {-# specialize fft'3 :: Array Int (Complex Float) -> Array Int (Complex Float) #-} {-# specialize fft'3 :: Array Int (Complex Double) -> Array Int (Complex Double) #-} -- | Length 3 FFT fft'3 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n] -> Array a (Complex b) -- ^ X[k] fft'3 a = array (0,2) [ (0, ((tmp1 + tmp4) :+ (tmp10 + tmp11))), (1, ((tmp5 + tmp8) :+ (tmp9 + tmp12))), (2, ((tmp5 - tmp8) :+ (tmp12 - tmp9))) ] where k866025403 = sqrt 3 / 2 k500000000 = 0.5 tmp1 = realPart (a!0) tmp10 = imagPart (a!0) tmp2 = realPart (a!1) tmp6 = imagPart (a!1) tmp3 = realPart (a!2) tmp7 = imagPart (a!2) tmp4 = tmp2 + tmp3 tmp9 = k866025403 * (tmp3 - tmp2) tmp8 = k866025403 * (tmp6 - tmp7) tmp11 = tmp6 + tmp7 tmp5 = tmp1 - (k500000000 * tmp4) tmp12 = tmp10 - (k500000000 * tmp11) {-# specialize fft'4 :: Array Int (Complex Float) -> Array Int (Complex Float) #-} {-# specialize fft'4 :: Array Int (Complex Double) -> Array Int (Complex Double) #-} -- | Length 4 FFT fft'4 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n] -> Array a (Complex b) -- ^ X[k] fft'4 a = array (0,3) [ (0, (tmp3 + tmp6) :+ (tmp15 + tmp16)), (1, (tmp11 + tmp14) :+ (tmp9 - tmp10)), (2, (tmp3 - tmp6) :+ (tmp15 - tmp16)), (3, (tmp11 - tmp14) :+ (tmp10 + tmp9)) ] where tmp1 = realPart (a!0) tmp7 = imagPart (a!0) tmp4 = realPart (a!1) tmp12 = imagPart (a!1) tmp2 = realPart (a!2) tmp8 = imagPart (a!2) tmp5 = realPart (a!3) tmp13 = imagPart (a!3) tmp3 = tmp1 + tmp2 tmp11 = tmp1 - tmp2 tmp9 = tmp7 - tmp8 tmp15 = tmp7 + tmp8 tmp6 = tmp4 + tmp5 tmp10 = tmp4 - tmp5 tmp14 = tmp12 - tmp13 tmp16 = tmp12 + tmp13 -------------------------------------------------------------------------------