-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.FIR.Kaiser
-- Copyright   :  (c) Matthew Donadio 1998
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module implements the Kaiser Window Method for designing FIR
-- filters.
--
-----------------------------------------------------------------------------

-- Reference:
--
-- @Book{dsp,
--   author = 	 "Alan V. Oppenheim and Ronald W. Schafer",
--   title = 	 "Discrete-Time Signal Processing",
--   publisher = 	 "Prentice-Hall",
--   year = 	 1989,
--   address =	 "Englewood Cliffs",
--   series =       {Prentice-Hall Signal Processing Series}
-- }

module DSP.Filter.FIR.Kaiser (kaiser_lpf, kaiser_hpf) where

import Data.Array

import DSP.Window
import DSP.Filter.FIR.Taps

-- Set the cutoff frequency to the middle of the transition band.  This
-- equation isn't numbered.

calc_wc :: Fractional a => a -> a -> a
calc_wc :: forall a. Fractional a => a -> a -> a
calc_wc a
wp a
ws = (a
wp forall a. Num a => a -> a -> a
+ a
ws) forall a. Fractional a => a -> a -> a
/ a
2

-- Equation 7.90

calc_dw :: Num a => a -> a -> a
calc_dw :: forall a. Num a => a -> a -> a
calc_dw a
wp a
ws = forall a. Num a => a -> a
abs (a
ws forall a. Num a => a -> a -> a
- a
wp)

-- Equation 7.91

calc_A :: (Floating a, Ord a) => a -> a -> a
calc_A :: forall a. (Floating a, Ord a) => a -> a -> a
calc_A a
d1 a
d2 = -a
20 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a -> a
logBase a
10 (forall a. Ord a => a -> a -> a
min a
d1 a
d2)

-- xEquation 7.92

calc_beta :: (Ord a, Floating a) => a -> a
calc_beta :: forall a. (Ord a, Floating a) => a -> a
calc_beta a
a | a
a forall a. Ord a => a -> a -> Bool
> a
50    = a
0.1102 forall a. Num a => a -> a -> a
* (a
a forall a. Num a => a -> a -> a
- a
8.7)
            | a
a forall a. Ord a => a -> a -> Bool
>= a
21   = a
0.5842 forall a. Num a => a -> a -> a
* ((a
aforall a. Num a => a -> a -> a
-a
21) forall a. Floating a => a -> a -> a
** a
0.4) forall a. Num a => a -> a -> a
+ a
0.07886 forall a. Num a => a -> a -> a
* (a
aforall a. Num a => a -> a -> a
-a
21)
            | Bool
otherwise = a
0.0

-- Equation 7.93

calc_M :: (Integral b, RealFrac a) => a -> a -> b
calc_M :: forall b a. (Integral b, RealFrac a) => a -> a -> b
calc_M a
a a
dw = forall a b. (RealFrac a, Integral b) => a -> b
ceiling ((a
a forall a. Num a => a -> a -> a
- a
8) forall a. Fractional a => a -> a -> a
/ (a
2.285 forall a. Num a => a -> a -> a
* a
dw))

-- Procedure on pg 455.  We should really check the peak approximation
-- error and then increase M if necessary.

-- | Designs a lowpass Kaiser filter

kaiser_lpf :: Double -- ^ wp
	   -> Double -- ^ ws
	   -> Double -- ^ dp
	   -> Double -- ^ ds
	   -> Array Int Double -- ^ h[n]

kaiser_lpf :: Double -> Double -> Double -> Double -> Array Int Double
kaiser_lpf Double
wp Double
ws Double
d1 Double
d2 = Array Int Double -> Array Int Double -> Array Int Double
window (Double -> Int -> Array Int Double
kaiser Double
beta Int
m) (forall a b.
(Ix a, Integral a, Enum b, Floating b, Eq b) =>
b -> a -> Array a b
lpf Double
wc Int
m)
    where wc :: Double
wc = forall a. Fractional a => a -> a -> a
calc_wc Double
wp Double
ws
          dw :: Double
dw = forall a. Num a => a -> a -> a
calc_dw Double
wp Double
ws
          a :: Double
a = forall a. (Floating a, Ord a) => a -> a -> a
calc_A Double
d1 Double
d2
          beta :: Double
beta = forall a. (Ord a, Floating a) => a -> a
calc_beta Double
a
          m :: Int
m = forall b a. (Integral b, RealFrac a) => a -> a -> b
calc_M Double
a Double
dw

-- The weird case for m below is because highpass (or bandstop) filters
-- should only be Type I.  Linear phase forces a null at w=pi for Type II
-- filters, which doesn't fit well with these kinds of filters.  Again,
-- we should really check the peak approximation error and then increase
-- M (by two) if necessary.

-- | Designs a highpass Kaiser filter

kaiser_hpf :: Double -- ^ wp
	   -> Double -- ^ ws
	   -> Double -- ^ dp
	   -> Double -- ^ ds
	   -> Array Int Double -- ^ h[n]

kaiser_hpf :: Double -> Double -> Double -> Double -> Array Int Double
kaiser_hpf Double
wp Double
ws Double
d1 Double
d2 = Array Int Double -> Array Int Double -> Array Int Double
window (Double -> Int -> Array Int Double
kaiser Double
beta Int
m) (forall a b.
(Ix a, Integral a, Enum b, Floating b, Eq b) =>
b -> a -> Array a b
hpf Double
wc Int
m)
    where wc :: Double
wc = forall a. Fractional a => a -> a -> a
calc_wc Double
wp Double
ws
          dw :: Double
dw = forall a. Num a => a -> a -> a
calc_dw Double
wp Double
ws
          a :: Double
a = forall a. (Floating a, Ord a) => a -> a -> a
calc_A Double
d1 Double
d2
          beta :: Double
beta = forall a. (Ord a, Floating a) => a -> a
calc_beta Double
a
          m :: Int
m = forall b. Integral b => b -> b
ceilingEven (forall b a. (Integral b, RealFrac a) => a -> a -> b
calc_M Double
a Double
dw)

ceilingEven :: Integral b => b -> b
ceilingEven :: forall b. Integral b => b -> b
ceilingEven b
x = b
x forall a. Num a => a -> a -> a
+ forall a. Integral a => a -> a -> a
mod (-b
x) b
2