diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2010, Judah Jacobson
+
+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 Judah Jacobson 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/Numeric/FFT/Vector/Base.hsc b/Numeric/FFT/Vector/Base.hsc
new file mode 100644
--- /dev/null
+++ b/Numeric/FFT/Vector/Base.hsc
@@ -0,0 +1,250 @@
+-- | A basic interface between Vectors and the fftw library.
+module Numeric.FFT.Vector.Base(
+            -- * Transforms
+            Transform(..),
+            planOfType,
+            PlanType(..),
+            plan,
+            run,
+            -- * Plans
+            Plan(),
+            planInputSize,
+            planOutputSize,
+            execute,
+            executeM,
+            -- * Unsafe C stuff
+            CFlags,
+            CPlan,
+            -- * Normalization helpers
+            Scalable(..),
+            modifyInput,
+            modifyOutput,
+            constMultOutput,
+            multC,
+            unsafeModify,
+            ) where
+
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as MS
+import Data.Vector.Generic as V hiding (forM_)
+import Data.Vector.Generic.Mutable as M
+import Data.List as L
+import Control.Monad.Primitive (RealWorld,PrimMonad(..),
+            unsafePrimToPrim, unsafePrimToIO)
+import Control.Monad(forM_)
+import Foreign (Storable(..), Ptr, unsafePerformIO, FunPtr,
+                ForeignPtr, withForeignPtr, newForeignPtr)
+import Foreign.C (CInt, CUInt)
+import Data.Bits ( (.&.) )
+import Data.Complex(Complex(..))
+import Foreign.Storable.Complex()
+
+
+
+#include <fftw3.h>
+
+---------------------
+-- Creating FFTW plans
+
+-- First, the Transform flags:
+data PlanType = Estimate | Measure | Patient | Exhaustive
+data Preservation = PreserveInput | DestroyInput
+
+type CFlags = CUInt
+
+-- | Marshal the Transform flags for use by fftw.
+planInitFlags :: PlanType -> Preservation -> CFlags
+planInitFlags pt pr = planTypeInt .&. preservationInt
+  where
+    planTypeInt = case pt of
+                    Estimate -> #const FFTW_ESTIMATE
+                    Measure -> #const FFTW_MEASURE
+                    Patient -> #const FFTW_PATIENT
+                    Exhaustive -> #const FFTW_EXHAUSTIVE
+    preservationInt = case pr of
+                    PreserveInput -> #const FFTW_PRESERVE_INPUT
+                    DestroyInput -> #const FFTW_DESTROY_INPUT
+
+newtype CPlan = CPlan {unCPlan :: ForeignPtr CPlan}
+
+withPlan :: CPlan -> (Ptr CPlan -> IO a) -> IO a
+withPlan = withForeignPtr . unCPlan
+
+foreign import ccall unsafe fftw_execute :: Ptr CPlan -> IO ()
+foreign import ccall "&" fftw_destroy_plan :: FunPtr (Ptr CPlan -> IO ())
+
+newPlan :: Ptr CPlan -> IO CPlan
+newPlan = fmap CPlan . newForeignPtr fftw_destroy_plan
+
+----------------------------------------
+-- vector-fftw plans
+
+-- | A 'Plan' can be used to run an @fftw@ algorithm for a specific input/output size.
+data Plan a b = Plan {
+                    planInput :: {-# UNPACK #-} !(VS.MVector RealWorld a),
+                    planOutput :: {-# UNPACK #-} !(VS.MVector RealWorld b),
+                    planExecute :: IO ()
+                }
+
+-- | The (only) valid input size for this plan.
+planInputSize :: Storable a => Plan a b -> Int
+planInputSize = MS.length . planInput
+
+-- | The (only) valid output size for this plan.
+planOutputSize :: Storable b => Plan a b -> Int
+planOutputSize = MS.length . planOutput
+
+-- | Run a plan on the given 'Vector'.
+--
+-- If @'planInputSize' p /= length v@, then calling
+-- @execute p v@ will throw an exception.
+execute :: (Vector v a, Vector v b, Storable a, Storable b) 
+            => Plan a b -> v a -> v b
+execute Plan{..} = \v -> -- fudge the arity to make sure it's always inlined
+    if n /= V.length v
+        then error $ "execute: size mismatch; expected " L.++ show n
+                    L.++ ", got " L.++ show (V.length v)
+        else unsafePerformIO $ do
+                        forM_ [0..n-1] $ \k -> M.unsafeWrite planInput k
+                                                $ V.unsafeIndex v k
+                        planExecute
+                        v' <- unsafeNew m
+                        forM_ [0..m-1] $ \k -> M.unsafeRead planOutput k
+                                                >>= M.unsafeWrite v' k
+                        V.unsafeFreeze v'
+  where
+    n = MS.length planInput
+    m = MS.length planOutput
+{-# INLINE execute #-}
+
+-- TODO: decide whether this is actually unsafe.
+-- | Run a plan on the given mutable vectors.  The same vector may be used for both
+-- input and output.
+--
+-- If @'planInputSize' p \/= length vIn@ or @'planOutputSize' p \/= length vOut@,
+-- then calling @unsafeExecuteM p vIn vOut@ will throw an exception.
+executeM :: forall m v a b . 
+        (PrimMonad m, MVector v a, MVector v b, Storable a, Storable b)
+            => Plan a b -- ^ The plan to run.
+            -> v (PrimState m) a  -- ^ The input vector.
+                    -> v (PrimState m) b -- ^ The output vector.
+                    -> m ()
+executeM Plan{..} = \vIn vOut ->
+    if n /= M.length vIn || m /= M.length vOut
+        then error $ "executeM: size mismatch; expected " L.++ show (n,m)
+                    L.++ ", got " L.++ show (M.length vIn, M.length vOut)
+        else unsafePrimToPrim $ act vIn vOut
+  where
+    n = MS.length planInput
+    m = MS.length planOutput
+
+    act :: v (PrimState m) a -> v (PrimState m) b -> IO ()
+    act vIn vOut = do
+            forM_ [0..n-1] $ \k -> unsafePrimToIO (M.unsafeRead vIn k :: m a)
+                                    >>= M.unsafeWrite planInput k
+            unsafePrimToPrim planExecute
+            forM_ [0..n-1] $ \k -> M.unsafeRead planOutput k
+                                    >>= unsafePrimToIO . (M.unsafeWrite vOut k
+                                                            :: b -> m ())
+{-# INLINE executeM #-}
+
+------------------
+-- Malloc/free of fftw array
+
+foreign import ccall unsafe fftw_malloc :: CInt -> IO (Ptr a)
+foreign import ccall "&" fftw_free :: FunPtr (Ptr a -> IO ())
+
+newFFTVector :: forall a . Storable a => Int -> IO (MS.MVector RealWorld a)
+newFFTVector n = do
+    p <- fftw_malloc $ toEnum $ n * sizeOf (undefined :: a)
+    fp <- newForeignPtr fftw_free p
+    return $ MS.MVector p n fp
+{-# INLINE newFFTVector #-}
+
+
+-----------------------
+-- Transforms: methods of plan creation.
+
+-- | A transform which may be applied to vectors of different sizes.
+data Transform a b = Transform {
+                        inputSize :: Int -> Int,
+                        outputSize :: Int -> Int,
+                        creationSizeFromInput :: Int -> Int,
+                        makePlan :: CInt -> Ptr a -> Ptr b -> CFlags -> IO (Ptr CPlan),
+                        normalization :: Int -> Plan a b -> Plan a b
+                    }
+
+-- | Create a 'Plan' of a specific size for this transform.
+planOfType :: (Storable a, Storable b) => PlanType
+                                -> Transform a b -> Int -> Plan a b
+planOfType ptype Transform{..} n
+  | m_in <= 0 || m_out <= 0 = error "Can't (yet) plan for empty arrays!"
+  | otherwise  = unsafePerformIO $ do
+    planInput <- newFFTVector m_in
+    planOutput <- newFFTVector m_out
+    MS.unsafeWith planInput $ \inP -> MS.unsafeWith planOutput $ \outP -> do
+    pPlan <- makePlan (toEnum n) inP outP $ planInitFlags ptype DestroyInput
+    cPlan <- newPlan pPlan
+    -- Use unsafeWith here to ensure that the Storable MVectors' ForeignPtrs
+    -- aren't released too soon:
+    let planExecute = MS.unsafeWith planInput $ \_ ->
+                        MS.unsafeWith planOutput $ \_ ->
+                          withPlan cPlan fftw_execute
+    return $ normalization n $ Plan {..}
+  where
+    m_in = inputSize n
+    m_out = outputSize n
+{-# INLINE planOfType #-}
+
+-- | Create a 'Plan' of a specific size.  This function is equivalent to
+-- @'planOfType' 'Estimate'@.
+plan :: (Storable a, Storable b) => Transform a b -> Int -> Plan a b
+plan = planOfType Estimate
+{-# INLINE plan #-}
+
+-- | Create and run a 'Plan' for the given transform.
+run :: (Vector v a, Vector v b, Storable a, Storable b)
+            => Transform a b -> v a -> v b
+run p = \v -> execute
+            (planOfType Estimate p $ creationSizeFromInput p $ V.length v)
+            v
+{-# INLINE run #-}
+
+---------------------------
+-- For scaling input/output:
+
+class Scalable a where
+    scaleByD :: Double -> a -> a
+    {-# INLINE scaleByD #-}
+
+instance Scalable Double where
+    scaleByD = (*)
+    {-# INLINE scaleByD #-}
+
+instance Scalable (Complex Double) where
+    scaleByD s (x:+y) = s*x :+ s*y
+    {-# INLINE scaleByD #-}
+
+
+{-# INLINE modifyInput #-}
+modifyInput :: (MS.MVector RealWorld a -> IO ()) -> Plan a b -> Plan a b
+modifyInput f p@Plan{..} = p {planExecute = f planInput >> planExecute}
+
+{-# INLINE modifyOutput #-}
+modifyOutput :: (MS.MVector RealWorld b -> IO ()) -> Plan a b -> Plan a b
+modifyOutput f p@Plan{..} = p {planExecute = planExecute >> f planOutput}
+
+{-# INLINE constMultOutput #-}
+constMultOutput :: (Storable b, Scalable b) => Double -> Plan a b -> Plan a b
+constMultOutput !s = modifyOutput (multC s)
+
+{-# INLINE multC #-}
+multC :: (Storable a, Scalable a) => Double -> MS.MVector RealWorld a -> IO ()
+multC !s v = forM_ [0..n-1] $ \k -> unsafeModify v k (scaleByD s)
+  where !n = MS.length v
+
+-- | Helper function; seems like it should be in the vector package...
+{-# INLINE unsafeModify #-}
+unsafeModify :: (Storable a)
+                => MS.MVector RealWorld a -> Int -> (a -> a) -> IO ()
+unsafeModify v k f = MS.unsafeRead v k >>= MS.unsafeWrite v k . f
diff --git a/Numeric/FFT/Vector/Invertible.hs b/Numeric/FFT/Vector/Invertible.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/FFT/Vector/Invertible.hs
@@ -0,0 +1,120 @@
+{- |
+This module provides normalized versions of the transforms in @fftw@.
+
+The forwards transforms in this module are identical to those in "Numeric.FFT.Vector.Unnormalized".
+The backwards transforms are normalized to be their inverse operations (approximately, due to floating point precision).
+
+For more information on the underlying transforms, see
+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.
+-}
+module Numeric.FFT.Vector.Invertible(
+                    -- * Creating and executing 'Plan's
+                    run,
+                    plan,
+                    execute,
+                    -- * Complex-to-complex transforms
+                    U.dft,
+                    idft,
+                    -- * Real-to-complex transforms
+                    U.dftR2C,
+                    dftC2R,
+                    -- * Real-to-real transforms
+                    -- $dct_size
+                    -- ** Discrete cosine transforms
+                    U.dct1,
+                    idct1,
+                    U.dct2,
+                    idct2,
+                    U.dct3,
+                    idct3,
+                    U.dct4,
+                    idct4,
+                    -- ** Discrete sine transforms
+                    U.dst1,
+                    idst1,
+                    U.dst2,
+                    idst2,
+                    U.dst3,
+                    idst3,
+                    U.dst4,
+                    idst4,
+                    ) where
+
+import Numeric.FFT.Vector.Base
+import qualified Numeric.FFT.Vector.Unnormalized as U
+import Data.Complex
+
+-- | A backward discrete Fourier transform which is the inverse of 'U.dft'.  The output and input sizes are the same (@n@).
+-- 
+-- @y_k = (1\/n) sum_(j=0)^(n-1) x_j e^(2pi i j k/n)@
+idft :: Transform (Complex Double) (Complex Double)
+idft = U.idft {normalization = \n -> constMultOutput $ 1 / toEnum n}
+
+-- | A normalized backward discrete Fourier transform which is the left inverse of
+-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)
+--
+-- This 'Transform' behaves differently than the others:
+--  
+--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
+--    /input/ size is @n \`div\` 2 + 1@.
+--
+--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
+--
+dftC2R :: Transform (Complex Double) Double
+dftC2R = U.dftC2R {normalization = \n -> constMultOutput $ 1 / toEnum n}
+
+-- OK, the inverse of each unnormalized operation.
+
+-- $dct_size
+-- The real-even (DCT) and real-odd (DST) transforms.  The input and output sizes
+-- are the same (@n@).
+
+
+-- | A type-1 discrete cosine transform which is the inverse of 'U.dct1'.
+--
+-- @y_k = (1\/(2(n-1)) [x_0 + (-1)^k x_(n-1) + 2 sum_(j=1)^(n-2) x_j cos(pi j k\/(n-1))]@
+idct1 :: Transform Double Double
+idct1 = U.dct1 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * (n-1))}
+
+-- | A type-3 discrete cosine transform which is the inverse of 'U.dct2'.
+-- 
+-- @y_k = (1\/(2n)) [x_0 + 2 sum_(j=1)^(n-1) x_j cos(pi j(k+1\/2)\/n)]@
+idct2 :: Transform Double Double
+idct2 = U.dct3 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
+-- | A type-2 discrete cosine transform which is the inverse of 'U.dct3'.
+-- 
+-- @y_k = (1\/n) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)k\/n)@
+idct3 :: Transform Double Double
+idct3 = U.dct2 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
+-- | A type-4 discrete cosine transform which is the inverse of 'U.dct4'.
+--
+-- @y_k = (1\/n) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)(k+1\/2)\/n)@
+idct4 :: Transform Double Double
+idct4 = U.dct4 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
+-- | A type-1 discrete sine transform which is the inverse of 'U.dst1'.
+--
+-- @y_k = (1\/(n+1)) sum_(j=0)^(n-1) x_j sin(pi(j+1)(k+1)\/(n+1))@
+idst1 :: Transform Double Double
+idst1 = U.dst1 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * (n+1))}
+
+-- | A type-3 discrete sine transform which is the inverse of 'U.dst2'.
+--
+-- @y_k = (1\/(2n)) [(-1)^k x_(n-1) + 2 sum_(j=0)^(n-2) x_j sin(pi(j+1)(k+1\/2)/n)]@
+idst2 :: Transform Double Double
+idst2 = U.dst3 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
+-- | A type-2 discrete sine transform which is the inverse of 'U.dst3'.
+--
+-- @y_k = (1\/n) sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1)\/n)@
+idst3 :: Transform Double Double
+idst3 = U.dst2 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
+-- | A type-4 discrete sine transform which is the inverse of 'U.dst4'.
+--
+-- @y_k = (1\/(2n)) sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1\/2)\/n)@
+idst4 :: Transform Double Double
+idst4 = U.dst4 {normalization = \n -> constMultOutput $ 1 / toEnum (2 * n)}
+
diff --git a/Numeric/FFT/Vector/Plan.hs b/Numeric/FFT/Vector/Plan.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/FFT/Vector/Plan.hs
@@ -0,0 +1,16 @@
+module Numeric.FFT.Vector.Plan(
+                -- * Transform
+                Transform(),
+                planOfType,
+                PlanType(..),
+                plan,
+                run,
+                -- * Plans
+                Plan(),
+                planInputSize,
+                planOutputSize,
+                execute,
+                executeM,
+                ) where
+
+import Numeric.FFT.Vector.Base
diff --git a/Numeric/FFT/Vector/Unitary.hs b/Numeric/FFT/Vector/Unitary.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/FFT/Vector/Unitary.hs
@@ -0,0 +1,126 @@
+{- |
+This module provides normalized versions of the transforms in @fftw@.
+
+All of the transforms are normalized so that
+
+ - Each transform is unitary, i.e., preserves the inner product and the sum-of-squares norm of its input.
+
+ - Each backwards transform is the inverse of the corresponding forwards transform.
+
+(Both conditions only hold approximately, due to floating point precision.)
+
+For more information on the underlying transforms, see
+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.
+-}
+module Numeric.FFT.Vector.Unitary(
+                -- * Creating and executing 'Plan's
+                run,
+                plan,
+                execute,
+                -- * Complex-to-complex transforms
+                dft,
+                idft,
+                -- * Real-to-complex transforms
+                dftR2C,
+                dftC2R,
+                -- * Discrete cosine transforms
+                -- $dct_size
+                dct2,
+                idct2,
+                dct4,
+                ) where
+
+import Numeric.FFT.Vector.Base
+import qualified Numeric.FFT.Vector.Unnormalized as U
+import Data.Complex
+import qualified Data.Vector.Storable.Mutable as MS
+import Control.Monad.Primitive(RealWorld)
+
+-- | A discrete Fourier transform. The output and input sizes are the same (@n@).
+--
+-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(-2pi i j k\/n)@
+dft :: Transform (Complex Double) (Complex Double)
+dft = U.dft {normalization = \n -> constMultOutput $ 1 / sqrt (toEnum n)}
+
+-- | An inverse discrete Fourier transform.  The output and input sizes are the same (@n@).
+-- 
+-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j e^(2pi i j k\/n)@
+idft :: Transform (Complex Double) (Complex Double)
+idft = U.idft {normalization = \n -> constMultOutput $ 1 / sqrt (toEnum n)}
+
+-- | A forward discrete Fourier transform with real data.  If the input size is @n@,
+-- the output size will be @n \`div\` 2 + 1@.
+dftR2C :: Transform Double (Complex Double)
+dftR2C = U.dftR2C {normalization = \n -> modifyOutput $
+                    complexR2CScaling (sqrt 2) n
+        }
+
+-- | A normalized backward discrete Fourier transform which is the left inverse of
+-- 'U.dftR2C'.  (Specifically, @run dftC2R . run dftR2C == id@.)
+--
+-- This 'Transform' behaves differently than the others:
+--  
+--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
+--    /input/ size is @n \`div\` 2 + 1@.
+--
+--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
+--
+dftC2R :: Transform (Complex Double) Double
+dftC2R = U.dftC2R {normalization = \n -> modifyInput $
+                    complexR2CScaling (sqrt 0.5) n
+        }
+
+complexR2CScaling :: Double -> Int -> MS.MVector RealWorld (Complex Double) -> IO ()
+complexR2CScaling !t !n !a = do
+    let !s1 = sqrt (1/toEnum n)
+    let !s2 = t * s1
+    let len = MS.length a
+    -- Justification for the use of unsafeModify:
+    -- The output size is 2n+1; so if n>0 then the output size is >=1;
+    -- and if n even then the output size is >=3.
+    unsafeModify a 0 $ scaleByD s1
+    if odd n
+        then multC s2 (MS.unsafeSlice 1 (len-1) a)
+        else do
+            unsafeModify a (len-1) $ scaleByD s1
+            multC s2 (MS.unsafeSlice 1 (len-2) a)
+
+
+-- $dct_size
+-- Some normalized real-even (DCT).  The input and output sizes
+-- are the same (@n@).
+
+
+-- | A type-4 discrete cosine transform.  It is its own inverse.
+-- 
+-- @y_k = (1\/sqrt n) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)(k+1\/2)\/n)@
+dct4 :: Transform Double Double
+dct4 = U.dct4 {normalization = \n -> constMultOutput $ 1 / sqrt (2 * toEnum n)}
+
+-- | A type-2 discrete cosine transform.  Its inverse is 'dct3'.
+--
+-- @y_k = w(k) sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)k\/n);@
+-- where
+-- @w(0)=1\/sqrt n@, and @w(k)=sqrt(2\/n)@ for @k>0@.
+dct2 :: Transform Double Double
+dct2 = U.dct2 {normalization = \n -> modifyOutput $ \a -> do
+    let n' = toEnum n
+    let !s1 = sqrt $ 1 / (4*n')
+    let !s2 = sqrt $ 1 / (2*n')
+    unsafeModify a 0 (*s1)
+    multC s2 (MS.unsafeSlice 1 (MS.length a-1) a)
+    }
+
+-- | A type-3 discrete cosine transform which is the inverse of 'dct2'.
+--
+-- @y_k = (-1)^k w(n-1) x_(n-1) + 2 sum_(j=0)^(n-2) w(j) x_j sin(pi(j+1)(k+1\/2)/n);@
+-- where
+-- @w(0)=1\/sqrt(n)@, and @w(k)=1/sqrt(2n)@ for @k>0@.
+idct2 :: Transform Double Double
+idct2 = U.dct3 {normalization = \n -> modifyInput $ \a -> do
+    let n' = toEnum n
+    let !s1 = sqrt $ 1 / n'
+    let !s2 = sqrt $ 1 / (2*n')
+    unsafeModify a 0 (*s1)
+    multC s2 (MS.unsafeSlice 1 (MS.length a-1) a)
+    }
diff --git a/Numeric/FFT/Vector/Unnormalized.hsc b/Numeric/FFT/Vector/Unnormalized.hsc
new file mode 100644
--- /dev/null
+++ b/Numeric/FFT/Vector/Unnormalized.hsc
@@ -0,0 +1,169 @@
+{- |
+Raw, unnormalized versions of the transforms in @fftw@.
+
+Note that the forwards and backwards transforms of this module are not actually
+inverses.  For example, @run idft (run dft v) /= v@ in general.
+
+For more information on the individual transforms, see
+<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.
+-}
+module Numeric.FFT.Vector.Unnormalized(
+                    -- * Creating and executing 'Plan's
+                    run,
+                    plan,
+                    execute,
+                    -- * Complex-to-complex transforms
+                    dft,
+                    idft,
+                    -- * Real-to-complex transforms
+                    dftR2C,
+                    dftC2R,
+                    -- * Real-to-real transforms
+                    -- $dct_size
+                    -- ** Discrete cosine transforms
+                    dct1,
+                    dct2,
+                    dct3,
+                    dct4,
+                    -- ** Discrete sine transforms
+                    dst1,
+                    dst2,
+                    dst3,
+                    dst4,
+                    ) where
+
+import Numeric.FFT.Vector.Base
+import Foreign
+import Foreign.C
+import Data.Complex
+
+#include <fftw3.h>
+
+-- | Whether the complex fft is forwards or backwards.
+type CDirection = CInt
+
+-- | The type of the cosine or sine transform.
+type CKind = (#type fftw_r2r_kind)
+
+foreign import ccall unsafe fftw_plan_dft_1d
+    :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> CDirection
+        -> CFlags -> IO (Ptr CPlan)
+
+foreign import ccall unsafe fftw_plan_dft_r2c_1d
+    :: CInt -> Ptr Double -> Ptr (Complex Double) -> CFlags -> IO (Ptr CPlan)
+
+foreign import ccall unsafe fftw_plan_dft_c2r_1d
+    :: CInt -> Ptr (Complex Double) -> Ptr Double -> CFlags -> IO (Ptr CPlan)
+
+foreign import ccall unsafe fftw_plan_r2r_1d
+    :: CInt -> Ptr Double -> Ptr Double -> CKind -> CFlags -> IO (Ptr CPlan)
+
+dft1D :: CDirection -> Transform (Complex Double) (Complex Double)
+dft1D d = Transform {
+            inputSize = id,
+            outputSize = id,
+            creationSizeFromInput = id,
+            makePlan = \n a b -> fftw_plan_dft_1d n a b d,
+            normalization = const id
+            }
+
+-- | A forward discrete Fourier transform.  The output and input sizes are the same (@n@).
+-- 
+-- @y_k = sum_(j=0)^(n-1) x_j e^(-2pi i j k/n)@
+dft :: Transform (Complex Double) (Complex Double)
+dft = dft1D (#const FFTW_FORWARD)
+
+-- | A backward discrete Fourier transform.  The output and input sizes are the same (@n@).
+-- 
+-- @y_k = sum_(j=0)^(n-1) x_j e^(2pi i j k/n)@
+idft :: Transform (Complex Double) (Complex Double)
+idft = dft1D (#const FFTW_BACKWARD)
+
+-- | A forward discrete Fourier transform with real data.  If the input size is @n@,
+-- the output size will be @n \`div\` 2 + 1@.
+dftR2C :: Transform Double (Complex Double)
+dftR2C = Transform {
+            inputSize = id,
+            outputSize = \n -> n `div` 2 + 1,
+            creationSizeFromInput = id,
+            makePlan = fftw_plan_dft_r2c_1d,
+            normalization = const id
+        }
+
+-- | A backward discrete Fourier transform which produces real data.
+--
+-- This 'Transform' behaves differently than the others:
+--  
+--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
+--    /input/ size is @n \`div\` 2 + 1@.
+--
+--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
+dftC2R :: Transform (Complex Double) Double
+dftC2R = Transform {
+            inputSize = \n -> n `div` 2 + 1,
+            outputSize = id,
+            creationSizeFromInput = \n -> 2 * (n-1),
+            makePlan = fftw_plan_dft_c2r_1d,
+            normalization = const id
+        }
+
+r2rTransform :: CKind -> Transform Double Double
+r2rTransform kind = Transform {
+                    inputSize = id,
+                    outputSize = id,
+                    creationSizeFromInput = id,
+                    makePlan = \n a b -> fftw_plan_r2r_1d n a b kind,
+                    normalization = const id
+                }
+
+-- $dct_size
+-- The real-even (DCT) and real-odd (DST) transforms.  The input and output sizes
+-- are the same (@n@).
+
+-- | A type-1 discrete cosine transform.  
+--
+-- @y_k = x_0 + (-1)^k x_(n-1) + 2 sum_(j=1)^(n-2) x_j cos(pi j k\/(n-1))@
+dct1 :: Transform Double Double
+dct1 = r2rTransform (#const  FFTW_REDFT00)
+
+-- | A type-2 discrete cosine transform.  
+--
+-- @y_k = 2 sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)k\/n)@
+dct2 :: Transform Double Double
+dct2 = r2rTransform (#const  FFTW_REDFT10)
+
+-- | A type-3 discrete cosine transform.  
+--
+-- @y_k = x_0 + 2 sum_(j=1)^(n-1) x_j cos(pi j(k+1\/2)\/n)@
+dct3 :: Transform Double Double
+dct3 = r2rTransform (#const  FFTW_REDFT01)
+
+-- | A type-4 discrete cosine transform.  
+--
+-- @y_k = 2 sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)(k+1\/2)\/n)@
+dct4 :: Transform Double Double
+dct4 = r2rTransform (#const  FFTW_REDFT11)
+
+-- | A type-1 discrete sine transform.
+-- 
+-- @y_k = 2 sum_(j=0)^(n-1) x_j sin(pi(j+1)(k+1)\/(n+1))@
+dst1 :: Transform Double Double
+dst1 = r2rTransform (#const  FFTW_RODFT00)
+
+-- | A type-2 discrete sine transform.
+-- 
+-- @y_k = 2 sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1)\/n)@
+dst2 :: Transform Double Double
+dst2 = r2rTransform (#const  FFTW_RODFT10)
+
+-- | A type-3 discrete sine transform.  
+--
+-- @y_k = (-1)^k x_(n-1) + 2 sum_(j=0)^(n-2) x_j sin(pi(j+1)(k+1\/2)/n)@
+dst3 :: Transform Double Double
+dst3 = r2rTransform (#const  FFTW_RODFT01)
+
+-- | A type-4 discrete sine transform.
+--
+-- @y_k = sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1\/2)\/n)@
+dst4 :: Transform Double Double
+dst4 = r2rTransform (#const FFTW_RODFT11)
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/vector-fftw.cabal b/vector-fftw.cabal
new file mode 100644
--- /dev/null
+++ b/vector-fftw.cabal
@@ -0,0 +1,46 @@
+Name:                vector-fftw
+
+Version:             0.1
+License:             BSD3
+License-file:        LICENSE
+Author:              Judah Jacobson
+Maintainer:          Judah Jacobson <judah.jacobson@gmail.com>
+Copyright:           (c) Judah Jacobson, 2010
+Category:            Math
+Build-type:          Simple
+Cabal-version:       >=1.6
+Synopsis:            A binding to the fftw library for one-dimensional vectors.
+Description:         This package provides bindings to the fftw library for one-dimensional vectors.
+                     It provides both high-level functions and more low-level manipulation of fftw plans.
+                     .
+                     We provide three different modules which wrap fftw's operations:
+                     .
+                      - "Numeric.FFT.Vector.Unnormalized" contains the raw transforms;
+                     . 
+                      - "Numeric.FFT.Vector.Invertible" scales the backwards transforms to be true inverses;
+                     .
+                      - "Numeric.FFT.Vector.Unitary" additionally scales all transforms to preserve the L2 (sum-of-squares) norm of the
+                        input.
+                     .
+                     Note that this package is currently not thread-safe.
+
+
+Library
+  Exposed-modules:     
+        Numeric.FFT.Vector.Unnormalized
+        Numeric.FFT.Vector.Invertible
+        Numeric.FFT.Vector.Unitary
+        Numeric.FFT.Vector.Plan
+
+  Other-modules:
+        Numeric.FFT.Vector.Base
+  
+  Build-depends: base==4.* && < 4.4, vector==0.7.*, primitive==0.3.*,
+                 storable-complex==0.2.*
+  Extra-libraries: fftw3
+
+  Extensions: ForeignFunctionInterface, RecordWildCards, BangPatterns, FlexibleInstances,
+                ScopedTypeVariables
+  ghc-options: -Wall
+
+  Ghc-Options: -O2