{- |
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>.
-}moduleNumeric.FFT.Vector.Unitary(-- * Creating and executing 'Plan'srun,plan,execute,-- * Complex-to-complex transformsdft,idft,-- * Real-to-complex transformsdftR2C,dftC2R,-- * Discrete cosine transforms-- $dct_sizedct2,idct2,dct4,)whereimportNumeric.FFT.Vector.BaseimportqualifiedNumeric.FFT.Vector.UnnormalizedasUimportData.CompleximportqualifiedData.Vector.Storable.MutableasMSimportControl.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(ComplexDouble)(ComplexDouble)dft=U.dft{normalization=\n->constMultOutput$1/sqrt(toEnumn)}-- | 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(ComplexDouble)(ComplexDouble)idft=U.idft{normalization=\n->constMultOutput$1/sqrt(toEnumn)}-- | A forward discrete Fourier transform with real data. If the input size is @n@,-- the output size will be @n \`div\` 2 + 1@.dftR2C::TransformDouble(ComplexDouble)dftR2C=U.dftR2C{normalization=\n->modifyOutput$complexR2CScaling(sqrt2)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(ComplexDouble)DoubledftC2R=U.dftC2R{normalization=\n->modifyInput$complexR2CScaling(sqrt0.5)n}complexR2CScaling::Double->Int->MS.MVectorRealWorld(ComplexDouble)->IO()complexR2CScaling!t!n!a=dolet!s1=sqrt(1/toEnumn)let!s2=t*s1letlen=MS.lengtha-- 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.unsafeModifya0$scaleByDs1ifoddnthenmultCs2(MS.unsafeSlice1(len-1)a)elsedounsafeModifya(len-1)$scaleByDs1multCs2(MS.unsafeSlice1(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::TransformDoubleDoubledct4=U.dct4{normalization=\n->constMultOutput$1/sqrt(2*toEnumn)}-- | 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::TransformDoubleDoubledct2=U.dct2{normalization=\n->modifyOutput$\a->doletn'=toEnumnlet!s1=sqrt$1/(4*n')let!s2=sqrt$1/(2*n')unsafeModifya0(*s1)multCs2(MS.unsafeSlice1(MS.lengtha-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::TransformDoubleDoubleidct2=U.dct3{normalization=\n->modifyInput$\a->doletn'=toEnumnlet!s1=sqrt$1/n'let!s2=sqrt$1/(2*n')unsafeModifya0(*s1)multCs2(MS.unsafeSlice1(MS.lengtha-1)a)}