From 31e73c0e56e57c7960487929f813016e8751f0ea Mon Sep 17 00:00:00 2001
From: Jesse McDonald <nybble41@gmail.com>
Date: Mon, 11 Nov 2013 04:29:54 -0600
Subject: [PATCH] Fix center value of kernel; refactor; add high-pass,
 band-reject, and band-pass filters.

---
 LowPass.hs | 100 +++++++++++++++++++++++++++++++++++------------------
 1 file changed, 67 insertions(+), 33 deletions(-)

diff --git a/LowPass.hs b/LowPass.hs
index a8acb0a..7ef88e1 100644
--- a/LowPass.hs
+++ b/LowPass.hs
@@ -12,32 +12,44 @@ import Numeric.GSL.Fourier          as LA
 
 import Control.Applicative
 import Data.Complex
+import System.IO
 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
+lowPassKernel :: Double -> Double -> Int -> Vector Double
+lowPassKernel sr fc ksize = raw / V.singleton sum
    where
-      n = V.enumFromN 0 kernel_size
-      t = (n - fromIntegral (kernel_size `div` 2)) / V.singleton sample_rate
+      n = V.enumFromN 0 ksize
+      t = (n - fromIntegral (ksize `div` 2)) / V.singleton sr
       -- 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)]
+      sinc' = sin (V.singleton (2*pi*fc) * t) / (V.singleton pi * t)
+      sinc  = sinc' V.// [(ksize `div` 2, 2 * fc)]
       -- Hamming window function
-      kmax = fromIntegral (kernel_size - 1)
+      kmax = fromIntegral (ksize - 1)
       hamm = 0.54 - 0.46 * cos (V.singleton (2 * pi / kmax) * n)
       -- Normalize the result
       raw = sinc * hamm
       sum = sumElements raw
 
+invertSpectrum :: Vector Double -> Vector Double
+invertSpectrum kernel = midVal `seq` (negate kernel V.// [(mid, 1 - midVal)])
+   where
+      mid = (dim kernel) `div` 2
+      midVal = kernel @> mid
+
+highPassKernel :: Double -> Double -> Int -> Vector Double
+highPassKernel sr fc ksize =
+   invertSpectrum $ lowPassKernel sr fc ksize
+
+bandRejectKernel :: Double -> (Double, Double) -> Int -> Vector Double
+bandRejectKernel sr (lfc, hfc) ksize =
+   lowPassKernel sr lfc ksize + highPassKernel sr hfc ksize
+
+bandPassKernel :: Double -> (Double, Double) -> Int -> Vector Double
+bandPassKernel sr (lfc, hfc) ksize =
+   invertSpectrum $ bandRejectKernel sr (lfc, hfc) ksize
+
 -- 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
@@ -48,18 +60,20 @@ convolve kernel input = V.generate osize $ \j -> rkernel <.> V.slice j ksize inp
       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
+decimate :: Int -> Vector Double -> Vector Double
+decimate osize vec =
+   V.generate osize $ \j ->
+      (sumElements $ V.slice (j * ssize) ssize vec) / fromIntegral ssize
+   where
+      vsize = dim vec
+      ssize = vsize `div` osize
+
+diffClockTimesSec :: ClockTime -> ClockTime -> Double
+diffClockTimesSec a b = sec + picosec / 1.0e12
+   where
+      diff    = diffClockTimes a b
+      sec     = fromIntegral $ tdSec diff
+      picosec = fromIntegral $ tdPicosec diff
 
 time :: IO a -> IO (a, Double)
 time f = do
@@ -67,10 +81,30 @@ time f = do
    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
+
+main :: IO ()
+main = do
+   let sample_rate = 10000 {-Hz-} :: Double
+   let cutoff      =  1000 {-Hz-} :: Double
+
+   let input_size  = 1000000 :: Int
+   let kernel_size =     201 :: Int
+
+   seed <- randomIO
+
+   (input,  inputTime)  <- time $ return $ LA.randomVector seed Gaussian input_size
+   (kernel, kernelTime) <- time $ return $ lowPassKernel sample_rate cutoff kernel_size
+   (result, resultTime) <- time $ return $ convolve kernel input
+
+   --V.mapM_ (printf "%10.6f\n") kernel
+
+   --let fft_result = V.map magnitude $ LA.fft $ V.map (:+0) result
+   --V.mapM_ (printf "%10.6f\n") . decimate 500 . V.take (dim fft_result `div` 2) $ fft_result
+
+   V.mapM_ (printf "%10.6f\n") . V.slice 0 50 $ result
+
+   hFlush stdout
+
+   hPutStrLn stderr $ printf "Input  Time: %8.6f seconds" inputTime
+   hPutStrLn stderr $ printf "Kernel Time: %8.6f seconds" kernelTime
+   hPutStrLn stderr $ printf "Result Time: %8.6f seconds" resultTime