-- | Discrete cosine transform (DCT-II).

dct :: U.Vector Double -> U.Vector Double

-dct = dct~~_~~ . G.map (:+0)

+dct = dctWorker . G.map (:+0)

--- | Discrete cosine transform, with complex coefficients (DCT-II).

+-- | Discrete cosine transform (DCT-II). Only real part of vector is

+-- transformed, imaginary part is ignored.

dct_ :: U.Vector CD -> U.Vector Double

-dct_ xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)

+dct_ = dctWorker . G.map (\(i :+ _) -> i :+ 0)

+dctWorker :: U.Vector CD -> U.Vector Double

+ = G.map realPart $ G.zipWith (*) weights (fft interleaved)

interleaved = G.backpermute xs $ G.enumFromThenTo 0 2 (len-2) G.++

G.enumFromThenTo (len-1) (len-3) 1

-- | Inverse discrete cosine transform (DCT-III). It's inverse of

-- 'dct' only up to scale parameter:

-- > (idct . dct) x = (* lenngth x)

idct :: U.Vector Double -> U.Vector Double

-idct = idct~~_~~ . G.map (:+0)

+idct = idctWorker . G.map (:+0)

--- | Inverse discrete cosine transform, with complex coefficients

+-- | Inverse discrete cosine transform (DCT-III). Only real part of vector is

+-- transformed, imaginary part is ignored.

idct_ :: U.Vector CD -> U.Vector Double

-idct_ xs = G.generate len interleave

+idct_ = idctWorker . G.map (\(i :+ _) -> i :+ 0)

+idctWorker :: U.Vector CD -> U.Vector Double

+idctWorker xs = G.generate len interleave

interleave z | even z = vals `G.unsafeIndex` halve z

| otherwise = vals `G.unsafeIndex` (len - halve z - 1)

-- | Inverse fast Fourier transform.

ifft :: U.Vector CD -> U.Vector CD

ifft xs = G.map ((/fi (G.length xs)) . conjugate) . fft . G.map conjugate $ xs