diff --git a/Data/Array/Repa/FFTW.hs b/Data/Array/Repa/FFTW.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/FFTW.hs
@@ -0,0 +1,180 @@
+{-|
+Module      : $Header$
+CopyRight   : (c) 8c6794b6, 2011, 2012
+License     : BSD3
+Maintainer  : 8c6794b6@gmail.com
+Stability   : experimental
+Portability : portable
+
+Performs fft of repa array data via FFTW.
+
+Currently supporting ('Complex' Double) arrays for dimensions 'DIM1', 'DIM2',
+and 'DIM3' only.
+
+-}
+
+module Data.Array.Repa.FFTW
+  ( -- * Examples
+    -- $examples
+
+    -- * Multi dimension functions
+    -- $multi
+
+    -- * References
+    -- $references
+
+    -- * FFT functions (1 dimension)
+    fft
+  , ifft
+    -- * FFT functions (2 dimension)
+  , fft2d
+  , ifft2d
+    -- * FFT functions (3 dimension)
+  , fft3d
+  , ifft3d
+  ) where
+
+import Data.Complex (Complex(..))
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Storable (Storable(..))
+import System.IO.Unsafe (unsafePerformIO)
+
+import Data.Array.CArray (CArray)
+import Data.Array.Repa ((:.)(..), Array, DIM1, DIM2, DIM3, Z(..))
+import Data.Array.Repa.Repr.ForeignPtr (F)
+import Foreign.Storable.Complex ()
+
+import qualified Data.Array.CArray as C
+import qualified Data.Array.Repa as R
+import qualified Data.Array.Repa.Repr.ForeignPtr as RF
+import qualified Math.FFT as FFT
+
+{-$examples
+
+Sample module:
+
+> import Data.Complex
+> import Data.Array.Repa
+> import Data.Array.Repa.Eval
+> import Data.Array.Repa.Repr.ForeignPtr
+> import Data.Array.Repa.FFTW
+>
+> a :: Array F DIM1 (Complex Double)
+> a = fromList (Z :. 4) [i :+ 0 | i <- [0..3]]
+
+Loading above in ghci:
+
+>>> toList a
+[0.0 :+ 0.0,1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]
+>>> toList $ fft a
+[6.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
+>>> toList $ ifft $ fft a
+[0.0 :+ 0.0,1.0 :+ 0.0,2.0 :+ 0.0,3.0 :+ 0.0]
+
+-}
+
+{-$multi
+
+Although FFTW library has more flexiblity, choice of dimensions in
+multi-dimension FFT functions ('fft2d', 'ifft2d' ... etc) were using fixed
+dimensions to perform FFTs. Those choices were made after @fft2dP@, @fft3dP@
+functions from @repa-algorithms@ package.
+
+-}
+
+{-$references
+
+* fftw : <http://www.fftw.org>
+
+* fftw haskell binding : <http://hackage.haskell.org/package/fft>
+
+* repa-algorithms: <http://hackage.haskell.org/package/repa-algorithms>
+
+-}
+
+-- --------------------------------------------------------------------------
+-- Exposed functions
+
+-- | Performs 1 dimension forward fft.
+fft :: Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double)
+fft = c2r . FFT.dft . r2c
+{-# INLINE fft #-}
+
+-- | Performs 1 dimension inverse fft.
+ifft :: Array F DIM1 (Complex Double) -> Array F DIM1 (Complex Double)
+ifft = c2r . FFT.idft . r2c
+{-# INLINE ifft #-}
+
+-- | Performs 2 dimension forward fft.
+fft2d :: Array F DIM2 (Complex Double) -> Array F DIM2 (Complex Double)
+fft2d = c2r2d . FFT.dftN [0,1] . r2c2d
+{-# INLINE fft2d #-}
+
+-- | Performs 2 dimension inverse fft.
+ifft2d :: Array F DIM2 (Complex Double) -> Array F DIM2 (Complex Double)
+ifft2d = c2r2d . FFT.idftN [0,1] . r2c2d
+{-# INLINE ifft2d #-}
+
+-- | Performs 3 dimension forward fft.
+fft3d :: Array F DIM3 (Complex Double) -> Array F DIM3 (Complex Double)
+fft3d = c2r3d . FFT.dftN [0,1,2] . r2c3d
+{-# INLINE fft3d #-}
+
+-- | Performs 3 dimension inverse fft.
+ifft3d :: Array F DIM3 (Complex Double) -> Array F DIM3 (Complex Double)
+ifft3d = c2r3d . FFT.idftN [0,1,2] . r2c3d
+{-# INLINE ifft3d #-}
+
+-- --------------------------------------------------------------------------
+-- Guts
+
+r2c :: Array F DIM1 (Complex Double) -> CArray Int (Complex Double)
+r2c rarr = unsafePerformIO $ do
+  let _:.nelem = R.extent rarr
+      fptr = RF.toForeignPtr rarr
+  C.unsafeForeignPtrToCArray fptr (0,nelem-1)
+{-# INLINE r2c #-}
+
+c2r :: CArray Int (Complex Double) -> Array F DIM1 (Complex Double)
+c2r carr = case C.toForeignPtr carr of
+  (n, fptr) -> let sh = Z:.n in
+    R.computeS $ R.fromFunction sh $ \ix ->
+    unsafePerformIO $ withForeignPtr fptr $ \ptr ->
+    peekElemOff ptr $ R.toIndex sh ix
+{-# INLINE c2r #-}
+
+r2c2d :: Array F DIM2 (Complex Double) -> CArray (Int, Int) (Complex Double)
+r2c2d rarr = unsafePerformIO $ do
+    let _:.n1:.n2 = R.extent rarr
+        fptr = RF.toForeignPtr rarr
+    C.unsafeForeignPtrToCArray fptr ((0,0), (n1-1, n2-1))
+{-# INLINE r2c2d #-}
+
+c2r2d :: CArray (Int, Int) (Complex Double) -> Array F DIM2 (Complex Double)
+c2r2d carr = case C.toForeignPtr carr of
+    (n, fptr) ->
+        let sh = Z:.n':.n'
+            n' = ceiling $ (sqrt $ fromIntegral n :: Double)
+        in  R.computeS $ R.fromFunction sh $ \ix ->
+            unsafePerformIO $ withForeignPtr fptr $ \ptr ->
+            peekElemOff ptr $ R.toIndex sh ix
+{-# INLINE c2r2d #-}
+
+r2c3d :: Array F DIM3 (Complex Double)
+      -> CArray (Int, Int, Int) (Complex Double)
+r2c3d rarr = unsafePerformIO $ do
+    let _:.n1:.n2:.n3 = R.extent rarr
+        fptr = RF.toForeignPtr rarr
+    C.unsafeForeignPtrToCArray fptr ((0,0,0), (n1-1, n2-1, n3-1))
+{-# INLINE r2c3d #-}
+
+c2r3d :: CArray (Int, Int, Int) (Complex Double)
+      -> Array F DIM3 (Complex Double)
+c2r3d carr = case C.toForeignPtr carr of
+    (n, fptr) ->
+        let sh = Z:.n':.n':.n'
+            n' = ceiling $ fromIntegral n ** (1/3 :: Double)
+        in  R.computeS $ R.fromFunction sh $ \ix ->
+            unsafePerformIO $ withForeignPtr fptr $ \ptr ->
+            peekElemOff ptr $ R.toIndex sh ix
+{-# INLINE c2r3d #-}
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2011-2013, 8c6794b6
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of 8c6794b6 nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/exec/bench.hs b/exec/bench.hs
new file mode 100644
--- /dev/null
+++ b/exec/bench.hs
@@ -0,0 +1,49 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-|
+Module      : $Header$
+CopyRight   : (c) 8c6794b6, 2011, 2012
+License     : BSD3
+Maintainer  : 8c6794b6@gmail.com
+Stability   : experimental
+Portability : portable
+
+Benchmark for comparing FFT with repa-fftw to repa-algorithms.
+-}
+module Main where
+
+import Control.DeepSeq (NFData(..))
+import Data.Complex
+import System.Random
+
+import Criterion.Main
+
+import Data.Array.Repa ((:.)(..), Array, U, DIM1, Z(..))
+import Data.Array.Repa.Repr.ForeignPtr (F)
+import qualified Data.Array.Repa as R
+import qualified Data.Array.Repa.Eval as Eval
+import qualified Data.Array.Repa.Algorithms.FFT as FFT
+import qualified Data.Array.Repa.FFTW as FFTW
+
+type RFFTArr = Array U DIM1 (Double, Double)
+
+main :: IO ()
+main = do
+  rs <- randomRs (0,1) `fmap` newStdGen
+  is <- randomRs (0,1) `fmap` newStdGen
+  let bench_fft n =
+        let mkarr ks = Eval.fromList (Z:.n) $ take n ks
+            ts = R.fromListUnboxed (Z:.n) $ take n $ zip rs is
+            cs = mkarr $ zipWith (:+) rs is
+            ra = FFT.fft1dP FFT.Forward :: RFFTArr -> IO RFFTArr
+        in  ts `seq` cs `seq`
+            bgroup ("n="++show n)
+              [ bench "repa-algorithms" (nfIO (ra ts))
+              , bench "repa-fftw" (nf FFTW.fft cs)
+              ]
+  defaultMain $ map (\k -> bench_fft (2^k)) [9..13::Int]
+
+instance NFData (Array U DIM1 (Double, Double)) where
+  rnf arr = arr `R.deepSeqArray` ()
+
+instance NFData (Array F DIM1 (Complex Double)) where
+  rnf arr = arr `R.deepSeqArray` ()
diff --git a/exec/test.hs b/exec/test.hs
new file mode 100644
--- /dev/null
+++ b/exec/test.hs
@@ -0,0 +1,94 @@
+-- | Main entry point for test executable.
+
+module Main where
+
+import Data.Array.Repa ((:.)(..), Array, DIM1, DIM2, DIM3, Shape, Z(..), toList)
+import Data.Array.Repa.Eval (fromList)
+import Data.Array.Repa.Repr.ForeignPtr (F)
+import Data.Complex (Complex(..), magnitude)
+
+import Test.Tasty (defaultMain, testGroup)
+import Test.Tasty.HUnit (Assertion, assertBool, testCase)
+import Test.Tasty.QuickCheck (Property, choose, forAll, testProperty, vector)
+
+import Data.Array.Repa.FFTW (fft, ifft, fft2d, ifft2d, fft3d, ifft3d)
+
+
+main :: IO ()
+main = defaultMain $
+    testGroup "repa-fftw"
+    [ testCase "magsum1d-hunit" case_magsum1d
+    , testCase "magsum2d-hunit" case_magsum2d
+    , testCase "magsum3d-hunit" case_magsum3d
+    , testProperty "magsum1d-qc" prop_magsum1d
+    , testProperty "magsum2d-qc" prop_magsum2d
+    , testProperty "magsum3d-qc" prop_magsum3d
+    ]
+
+case_magsum1d :: Assertion
+case_magsum1d = do
+    let a = fromList (Z :. 4) [i :+ 0 | i <- [0..3]]
+    assertBool "diff is greater than epsilon" $ compareMagnitudeSum1d a
+
+case_magsum2d :: Assertion
+case_magsum2d = do
+    let a = fromList (Z :. 4 :. 4) [i :+ 0 | i <- [0..15]]
+    assertBool "diff is greater than epsilon" $ compareMagnitudeSum2d a
+
+case_magsum3d :: Assertion
+case_magsum3d = do
+    let a = fromList (Z :. 4 :. 4 :. 4) [i :+ 0 | i <- [0..63]]
+    assertBool "diff is greater than epsilon" $ compareMagnitudeSum3d a
+
+prop_magsum1d :: Property
+prop_magsum1d =
+    forAll (choose (1, 16::Int)) $ \pow ->
+    let size = (2 ^ pow)
+    in forAll (vector size) $ \vals ->
+       compareMagnitudeSum1d $ fromList (Z :. size) vals
+
+prop_magsum2d :: Property
+prop_magsum2d =
+    forAll (choose (1, 4::Int)) $ \pow ->
+    let extnt = 2 ^ pow
+        size  = extnt ^ (2 :: Int)
+    in  forAll (vector size) $ \vals ->
+        compareMagnitudeSum2d $ fromList (Z :. extnt :. extnt) vals
+
+prop_magsum3d :: Property
+prop_magsum3d =
+    forAll (choose (1, 4::Int)) $ \pow ->
+    let extnt = 2 ^ pow
+        size  = extnt ^ (3 :: Int)
+    in  forAll (vector size) $ \vals ->
+        compareMagnitudeSum3d $ fromList (Z :. extnt :. extnt :. extnt) vals
+
+-- | Performs fft and ifft, then compareing magnitude of sum with the
+-- original.
+compareMagnitudeSum1d :: Array F DIM1 (Complex Double) -> Bool
+compareMagnitudeSum1d = compareMagnitudeSum fft ifft
+
+-- | Like 'compareMagnitudeSum', for 2-dimensional array.
+compareMagnitudeSum2d :: Array F DIM2 (Complex Double) -> Bool
+compareMagnitudeSum2d = compareMagnitudeSum fft2d ifft2d
+
+-- | Like 'compareMagnitudeSum', for 3-dimensional array.
+compareMagnitudeSum3d :: Array F DIM3 (Complex Double) -> Bool
+compareMagnitudeSum3d = compareMagnitudeSum fft3d ifft3d
+
+-- | Apply given functions, then compare the result with original.
+compareMagnitudeSum ::
+    (Shape sh)
+    => (Array F sh (Complex Double) -> a)
+    -- ^ Forward FFT function.
+    -> (a -> Array F sh (Complex Double))
+    -- ^ Inverse FFT function.
+    -> Array F sh (Complex Double)
+    -- ^ Original input.
+    -> Bool
+compareMagnitudeSum fn ifn arr = diff < epsilon
+  where
+    diff = abs (mag arr - mag arr')
+    mag = magnitude . sum . toList
+    arr' = ifn $ fn arr
+    epsilon = 1e-2
diff --git a/repa-fftw.cabal b/repa-fftw.cabal
new file mode 100644
--- /dev/null
+++ b/repa-fftw.cabal
@@ -0,0 +1,72 @@
+Name:                repa-fftw
+Version:             3.2.3.1
+Synopsis:            Perform fft with repa via FFTW
+License:             BSD3
+License-file:        LICENSE
+Author:              <8c6794b6@gmail.com>
+Maintainer:          <8c6794b6@gmail.com>
+Category:            Math
+Build-type:          Simple
+Cabal-version:       >=1.8
+
+Description:
+  Performs FFT on repa array data with fftw.
+
+flag dev
+  description:
+    Activate developer flags
+  default:
+    False
+
+source-repository head
+  type: git
+  location: https://github.com/8c6794b6/repa-fftw.git
+
+Library
+  ghc-options:
+    -Wall -fno-warn-orphans
+    -O3 -fllvm -optl-O3
+  exposed-modules:
+    Data.Array.Repa.FFTW
+  build-depends:
+    base >= 4.6 && < 5,
+    carray >= 0.1.5 && < 0.2,
+    fft >= 0.1.5 && < 0.2,
+    repa >= 3.2.3,
+    storable-complex >= 0.2
+
+test-suite test
+  type:
+    exitcode-stdio-1.0
+  hs-source-dirs:
+    exec
+  main-is:
+    test.hs
+  ghc-options:
+    -Wall
+  build-depends:
+    base >= 4.6 && < 5,
+    repa >= 3.2.3,
+    repa-fftw -any,
+    tasty >= 0.4.0,
+    tasty-hunit >= 0.4.1,
+    tasty-quickcheck >= 0.3.1
+
+benchmark bench
+  type:
+    exitcode-stdio-1.0
+  hs-source-dirs:
+    exec
+  main-is:
+    bench.hs
+  ghc-options:
+    -Wall -fno-warn-orphans
+    -threaded -rtsopts
+  build-depends:
+    base >= 4.6 && < 5,
+    criterion >= 0.8,
+    deepseq >= 1.3.0 && < 2,
+    repa >= 3.2.3,
+    repa-algorithms >= 3.2,
+    repa-fftw -any,
+    random >= 1.0.1.1