From 9a1803ec718b5a35dbbb5da00bfe6c1db53f897d Mon Sep 17 00:00:00 2001
From: Jesse McDonald <nybble41@gmail.com>
Date: Sun, 10 Nov 2013 05:41:23 -0600
Subject: [PATCH] Initial commit.

---
 .gitignore |  8 ++++++
 LowPass.hs | 72 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Makefile   | 19 ++++++++++++++
 3 files changed, 99 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 LowPass.hs
 create mode 100644 Makefile

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..edc097e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,8 @@
+LowPass
+dist
+cabal-dev
+*.o
+*.hi
+*.chi
+*.chs.h
+.virthualenv
diff --git a/LowPass.hs b/LowPass.hs
new file mode 100644
index 0000000..7ac0c0f
--- /dev/null
+++ b/LowPass.hs
@@ -0,0 +1,72 @@
+{-# 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 Control.Applicative
+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
+   (input,  inputTime)  <- time $ V.fromListN isize . randoms <$> newStdGen
+   (_,      kernelTime) <- time $ return lowPassKernel
+   (result, resultTime) <- time $ return $ convolve lowPassKernel input
+   V.mapM_ (printf "%0.6f\n") $ V.slice 0 50 result
+   printf "Input  Time: %0.6f seconds\n" $ inputTime
+   printf "Kernel Time: %0.6f seconds\n" $ kernelTime
+   printf "Result Time: %0.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
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..344c1cf
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,19 @@
+GHC_OPTS = -Odph -rtsopts -threaded -fno-liberate-case -funfolding-use-threshold1000 \
+           -funfolding-keeness-factor1000 -fllvm -optlo-O3 \
+	   -package hmatrix -package old-time
+
+all: LowPass
+
+clean:
+	$(RM) LowPass LowPass.o LowPass.hi
+
+LowPass: LowPass.o
+	ghc $(GHC_OPTS) $^ -o $@
+
+%.o: %.hs
+	ghc $(GHC_OPTS) -c $< -o $@
+
+# DO NOT DELETE: Beginning of Haskell dependencies
+LowPass.o : LowPass.hs
+# DO NOT DELETE: End of Haskell dependencies
+