vector-fftw (empty) → 0.1
raw patch · 8 files changed
+759/−0 lines, 8 filesdep +basedep +primitivedep +storable-complexsetup-changed
Dependencies added: base, primitive, storable-complex, vector
Files
- LICENSE +30/−0
- Numeric/FFT/Vector/Base.hsc +250/−0
- Numeric/FFT/Vector/Invertible.hs +120/−0
- Numeric/FFT/Vector/Plan.hs +16/−0
- Numeric/FFT/Vector/Unitary.hs +126/−0
- Numeric/FFT/Vector/Unnormalized.hsc +169/−0
- Setup.hs +2/−0
- vector-fftw.cabal +46/−0
+ LICENSE view
@@ -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.
+ Numeric/FFT/Vector/Base.hsc view
@@ -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
+ Numeric/FFT/Vector/Invertible.hs view
@@ -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)}+
+ Numeric/FFT/Vector/Plan.hs view
@@ -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
+ Numeric/FFT/Vector/Unitary.hs view
@@ -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)+ }
+ Numeric/FFT/Vector/Unnormalized.hsc view
@@ -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)
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ vector-fftw.cabal view
@@ -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