Forward DCT implemented in Haskell. Familiarize myself how to use mutable vector in ST monad. DCT algorithm might contain bugs :). Modified from ijgjpeg fdctflt.

DctQuant.hs

module DctQuant where

import Prelude hiding (mapM_, read)
import Data.Vector.Mutable
import Data.Vector
import Control.Monad.ST
import Control.Monad.Primitive (PrimMonad, PrimState)

fdct :: PrimMonad m => Vector Float -> m (MVector (PrimState m) Float)
fdct vi = do
    v <- new 64
    flip mapM_ (fromList [0..7])
                  (\i ->
                  let elemptr = Data.Vector.slice i 8 vi
                      tmp0 = elemptr ! 0 + elemptr ! 7
                      tmp7 = elemptr ! 0 - elemptr ! 7
                      tmp1 = elemptr ! 1 + elemptr ! 6
                      tmp6 = elemptr ! 1 - elemptr ! 6
                      tmp2 = elemptr ! 2 + elemptr ! 5
                      tmp5 = elemptr ! 2 - elemptr ! 5
                      tmp3 = elemptr ! 3 + elemptr ! 4
                      tmp4 = elemptr ! 3 - elemptr ! 4

                      -- even part
                      tmp10 = tmp0 + tmp3
                      tmp13 = tmp0 - tmp3
                      tmp11 = tmp1 + tmp2
                      tmp12 = tmp1 - tmp2

                      z1 = (tmp12 + tmp13) * 0.707106781

                      -- odd part
                      tmp10' = tmp4 + tmp5
                      tmp11' = tmp5 + tmp6
                      tmp12' = tmp6 + tmp7

                      z5 = (tmp10'-tmp12') * 0.382683433
                      z2 = 0.541196100 * tmp10'+z5
                      z4 = 1.306562965 * tmp12'+z5
                      z3 = tmp11' * 0.707106781 :: Float

                      z11 = tmp7 + z3
                      z13 = tmp7 - z3 in
                      do
                        write v (i*8+0) (tmp10+tmp11-8.0*128) -- phase 3
                        write v (i*8+4) (tmp10-tmp11)
                        write v (i*8+2) (tmp13+z1)
                        write v (i*8+6) (tmp13-z1)
                        write v (i*8+5) (z13+z2)
                        write v (i*8+3) (z13-z2)
                        write v (i*8+1) (z11+z4)
                        write v (i*8+7) (z11-z4))

    -- phase 2: process columns
    flip mapM_ (fromList [0..7])
                  (\i ->
                   do
                     t0 <- read v (8*0+i)
                     t1 <- read v (8*1+i)
                     t2 <- read v (8*2+i)
                     t3 <- read v (8*3+i)
                     t4 <- read v (8*4+i)
                     t5 <- read v (8*5+i)
                     t6 <- read v (8*6+i)
                     t7 <- read v (8*7+i)
                     let tmp0 = t0+t7
                         tmp7 = t0-t7
                         tmp1 = t1+t6
                         tmp6 = t1-t6
                         tmp2 = t2+t5
                         tmp5 = t2-t5
                         tmp3 = t3+t4
                         tmp4 = t3-t4
                         -- even part
                         tmp10 = tmp0+tmp3
                         tmp13 = tmp0-tmp3
                         tmp11 = tmp1+tmp2
                         tmp12 = tmp1-tmp2
                         z1 = (tmp12+tmp13)*0.707106781 -- c4

                         tmp10' = tmp4+tmp5
                         tmp11' = tmp5+tmp6
                         tmp12' = tmp6+tmp7
                         z5 = (tmp10'-tmp12')*0.382683433
                         z2 = 0.541196100*tmp10'+z5
                         z4 = 1.306562965*tmp12'+z5
                         z3 = tmp11'*0.707106781
                         z11 = tmp7+z3
                         z13 = tmp7-z3

                     write v (8*0+i) (tmp10+tmp11)
                     write v (8*4+i) (tmp10-tmp11)
                     write v (8*2+i) (tmp13+z1)
                     write v (8*6+i) (tmp13-z1)
                     write v (8*5+i) (z13+z2)
                     write v (8*3+i) (z13-z2)
                     write v (8*1+i) (z11+z4)
                     write v (8*7+i) (z11-z4))
    return v

-- How to use:
-- let mat = Prelude.map (fromIntegral::Int->Float) [1,2,3,4,5,6,7,8,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16,
--                                                   9,10,11,12,13,14,15,16]

-- main = runST $ ((fdct $ fromList mat) >>= freeze)