-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.FIR.Smooth
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Herrmann type smooth FIR filters, from Hamming, Chapter 7, also
-- known as maximally flat FIR filters
--
-- If x is the -3 dB point, then p\/q = -(x+1)\/(x-1)
--
-----------------------------------------------------------------------------

-- TODO: function for rational fraction approximation

-- TODO: input parameters in the style of sect53.f

module DSP.Filter.FIR.Smooth (smoothfir) where

import Data.Array

import Polynomial.Basic

-- Normalize is the step to set g(1) = 1 (pg 123)

normalize :: Fractional a => [a] -> [a]
normalize :: forall a. Fractional a => [a] -> [a]
normalize [a]
x = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/ a
a) [a]
x
    where a :: a
a = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
x

-- Expand performs the algorithm in Sect 7.3

expand :: Fractional a => [a] -> [a]
expand :: forall a. Fractional a => [a] -> [a]
expand (a
x1:a
x2:[]) = [ a
x1, a
x2 ]
expand (a
x:[a]
xs) = forall a. Fractional a => a -> [a] -> [a]
expand' a
x forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => [a] -> [a]
expand [a]
xs

expand' :: Fractional a => a -> [a] -> [a]
expand' :: forall a. Fractional a => a -> [a] -> [a]
expand' a
x [a]
ys0 = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+) (a
x forall a. a -> [a] -> [a]
: forall a. Fractional a => [a] -> [a]
m1 [a]
ys0) (forall a. Fractional a => [a] -> [a]
p1 [a]
ys0)
    where m1 :: [a] -> [a]
m1 (a
y:[a]
ys) = a
y forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map (a
0.5forall a. Num a => a -> a -> a
*) [a]
ys
	  p1 :: [a] -> [a]
p1 (a
_:[a]
ys) = forall a b. (a -> b) -> [a] -> [b]
map (a
0.5forall a. Num a => a -> a -> a
*) [a]
ys forall a. [a] -> [a] -> [a]
++ [ a
0, a
0 ]

-- Reflect makes the filter symmetric (not sure where this is stated)

reflect :: Fractional a => [a] -> [a]
reflect :: forall a. Fractional a => [a] -> [a]
reflect (a
x:[a]
xs) = (forall a b. (a -> b) -> [a] -> [b]
map (a
0.5forall a. Num a => a -> a -> a
*) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse [a]
xs) forall a. [a] -> [a] -> [a]
++ a
x forall a. a -> [a] -> [a]
: (forall a b. (a -> b) -> [a] -> [b]
map (a
0.5forall a. Num a => a -> a -> a
*) [a]
xs)

-- The actual function.  Note that we use (1+t)^p * (1-t)^q directly
-- since we have a polynomial library.

-- | designs smooth FIR filters

smoothfir :: (Ix a, Integral a, Fractional b) => a -- ^ p
	  -> a -- ^ q
	  -> Array a b -- ^ h[n]

smoothfir :: forall a b. (Ix a, Integral a, Fractional b) => a -> a -> Array a b
smoothfir a
p a
q = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a
0,a
nforall a. Num a => a -> a -> a
-a
1) forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => [a] -> [a]
reflect forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => [a] -> [a]
expand forall a b. (a -> b) -> a -> b
$ [b]
b
    where b' :: [b]
b' = forall a. Num a => [a] -> [a] -> [a]
polymult (forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [ b
1, b
1 ] a
p) (forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [ b
1, -b
1 ] a
q)
          b1 :: [b]
b1 = forall a. Fractional a => [a] -> a -> [a]
polyinteg [b]
b' b
0
	  c :: b
c = -forall a. Num a => [a] -> a -> a
polyeval [b]
b1 (-b
1)
	  b :: [b]
b = forall a. Fractional a => [a] -> [a]
normalize forall a b. (a -> b) -> a -> b
$ b
c forall a. a -> [a] -> [a]
: forall a. [a] -> [a]
tail [b]
b1
	  n :: a
n = a
2 forall a. Num a => a -> a -> a
* (a
pforall a. Num a => a -> a -> a
+a
1 forall a. Num a => a -> a -> a
+ a
qforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
- a
1

-- Test

-- map (256*) $ elems $ smoothfir 3 1 == [ -1, -5, -5, 20, 70, 98, 70, 20, -5, -5, -1 ]