{-# LANGUAGE BangPatterns #-}

module Main (main) where

import Prelude                      as P
import Control.Monad                as M
import Data.List                    as L
import Data.Vector.Storable         as V
import Data.Vector.Storable.Mutable as MV
import Numeric.LinearAlgebra        as LA
import Numeric.GSL.Fourier          as LA

import Control.Applicative
import Data.Complex
import System.Random
import System.Time
import Text.Printf

sample_rate, cutoff :: Double
sample_rate = 10000 {-Hz-}
cutoff      =  1000 {-Hz-}

kernel_size :: Int
kernel_size = 101

lowPassKernel :: Vector Double
lowPassKernel = raw / V.singleton sum
   where
      n = V.enumFromN 0 kernel_size
      t = (n - fromIntegral (kernel_size `div` 2)) / V.singleton sample_rate
      -- sinc function; replace division by zero with limit when t=0
      sinc' = sin (V.singleton (2*pi*cutoff) * t) / (V.singleton pi * t)
      sinc  = sinc' V.// [(kernel_size `div` 2, 2 * pi * cutoff / sample_rate)]
      -- Hamming window function
      kmax = fromIntegral (kernel_size - 1)
      hamm = 0.54 - 0.46 * cos (V.singleton (2 * pi / kmax) * n)
      -- Normalize the result
      raw = sinc * hamm
      sum = sumElements raw

-- convolution = integral(kernel(t-tau)*input(tau),tau)
-- t is output vector index (j); products and summation are done with dot product (<.>)
convolve :: Vector Double -> Vector Double -> Vector Double
convolve kernel input = V.generate osize $ \j -> rkernel <.> V.slice j ksize input
   where
      ksize = dim rkernel
      isize = dim input
      osize = isize - ksize
      rkernel = V.reverse kernel

main :: IO ()
main = do
   let isize = 2000000
   seed <- randomIO
   (input,  inputTime)  <- time $ return $ LA.randomVector seed Gaussian isize
   (_,      kernelTime) <- time $ return lowPassKernel
   (result, resultTime) <- time $ return $ convolve lowPassKernel input
   V.mapM_ (printf "%10.6f\n") $ V.slice 0 50 result
   --V.mapM_ (printf "%10.6f\n") $ V.map magnitude $ LA.fft $ V.map (:+0) result
   printf "Input  Time: %8.6f seconds\n" $ inputTime
   printf "Kernel Time: %8.6f seconds\n" $ kernelTime
   printf "Result Time: %8.6f seconds\n" $ resultTime

time :: IO a -> IO (a, Double)
time f = do
   start <- getClockTime
   x     <- f
   end   <- x `seq` getClockTime
   return (x, diffClockTimesSec end start)
 
diffClockTimesSec :: ClockTime -> ClockTime -> Double
diffClockTimesSec a b = sec + picosec / 1.0e12
   where
      diff    = diffClockTimes a b
      sec     = fromIntegral $ tdSec diff
      picosec = fromIntegral $ tdPicosec diff