-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Window
-- Copyright   :  (c) Matthew Donadio 1998
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Commonly used window functions.  Except for the Parzen window, the
-- results of all of these /look/ right, but I have to check them against
-- either Matlab or my C code.
--
-- More windowing functions exist, but I have to dig through my papers to
-- find the equations.
--
-----------------------------------------------------------------------------

-- TODO: These functions should probably be reworked to use list
-- comprehensions...

{-

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}
}

@Book{kay,
  author =       "Steven M. Kay",
  title =        "Modern Spectral Estimation: Theory \& Application",
  publisher =    "Prentice Hall",
  year =         1988,
  address =      "Englewood Cliffs",
  series =       {Prentice-Hall Signal Processing Series}
}

-}

module DSP.Window
    (window, rectangular, bartlett, hanning, hamming, blackman,
     kaiser, gen_hamming, parzen) where

import DSP.Basic ((^!))
import Data.Array


-- | Applys a window, @w@, to a sequence @x@

window :: Array Int Double -- ^ w[n]
       -> Array Int Double -- ^ x[n]
       -> Array Int Double -- ^ w[n] * x[n]

window :: Array Int Double -> Array Int Double -> Array Int Double
window Array Int Double
w Array Int Double
x = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
m) [ Array Int Double
wforall i e. Ix i => Array i e -> i -> e
!Int
i forall a. Num a => a -> a -> a
* Array Int Double
xforall i e. Ix i => Array i e -> i -> e
!Int
i | Int
i <- [Int
0..Int
m] ]
    where m :: Int
m = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int Double
w

-- | rectangular window

rectangular :: Int -- ^ M
            -> Array Int Double -- ^ w[n]

rectangular :: Int -> Array Int Double
rectangular Int
m = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
m) forall a b. (a -> b) -> a -> b
$ forall a. a -> [a]
repeat Double
1

-- | Bartlett  window

bartlett :: Int -- ^ M
         -> Array Int Double -- ^ w[n]

bartlett :: Int -> Array Int Double
bartlett = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
bartlett'

bartlett' :: Double -> Double -> Double
bartlett' :: Double -> Double -> Double
bartlett' Double
m Double
n =
   if Double
n forall a. Ord a => a -> a -> Bool
<= Double
m forall a. Fractional a => a -> a -> a
/ Double
2
     then Double
2 forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m
     else Double
2 forall a. Num a => a -> a -> a
- Double
2 forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m

-- | Hanning window

hanning :: Int -- ^ M
        -> Array Int Double -- ^ w[n]

hanning :: Int -> Array Int Double
hanning = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
hanning'

hanning' :: Double -> Double -> Double
hanning' :: Double -> Double -> Double
hanning' Double
m Double
n = Double
0.5 forall a. Num a => a -> a -> a
- Double
0.5 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m)

-- | Hamming window

hamming :: Int -- ^ M
        -> Array Int Double -- ^ w[n]

hamming :: Int -> Array Int Double
hamming = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
hamming'

hamming' :: Double -> Double -> Double
hamming' :: Double -> Double -> Double
hamming' Double
m Double
n = Double
0.54 forall a. Num a => a -> a -> a
- Double
0.46 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m)


-- | Blackman window

blackman :: Int -- ^ M
         -> Array Int Double -- ^ w[n]

blackman :: Int -> Array Int Double
blackman = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
blackman'

blackman' :: Double -> Double -> Double
blackman' :: Double -> Double -> Double
blackman' Double
m Double
n =
   Double
0.42 forall a. Num a => a -> a -> a
- Double
0.5 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m) forall a. Num a => a -> a -> a
+
   Double
0.08 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (Double
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m)

-- | Generalized Hamming window

gen_hamming :: Double -- ^ alpha
            -> Int -- ^ M
            -> Array Int Double -- ^ w[n]

gen_hamming :: Double -> Int -> Array Int Double
gen_hamming = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double -> Double
gen_hamming'

gen_hamming' :: Double -> Double -> Double -> Double
gen_hamming' :: Double -> Double -> Double -> Double
gen_hamming' Double
a Double
m Double
n = Double
a forall a. Num a => a -> a -> a
- (Double
1 forall a. Num a => a -> a -> a
- Double
a) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
n forall a. Fractional a => a -> a -> a
/ Double
m)

-- | rectangular window

kaiser :: Double -- ^ beta
       -> Int -- ^ M
       -> Array Int Double -- ^ w[n]

kaiser :: Double -> Int -> Array Int Double
kaiser = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double -> Double -> Double
kaiser'

kaiser' :: Double -> Double -> Double -> Double
kaiser' :: Double -> Double -> Double -> Double
kaiser' Double
b Double
m Double
n =
   let a :: Double
a = Double
m forall a. Fractional a => a -> a -> a
/ Double
2
   in  Double -> Double
i0 (Double
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (Double
1 forall a. Num a => a -> a -> a
-((Double
nforall a. Num a => a -> a -> a
-Double
a)forall a. Fractional a => a -> a -> a
/Double
a)forall a. Num a => a -> Int -> a
^!Int
2)) forall a. Fractional a => a -> a -> a
/ Double -> Double
i0 Double
b


-- Recursive computation of I0, the zeroth-order modified Bessel function
-- of the first kind.

i0  :: Double -> Double
i0 :: Double -> Double
i0 Double
x = Double -> Double -> Double -> Double
i0' Double
x Double
2 Double
1

i0' :: Double -> Double -> Double -> Double
i0' :: Double -> Double -> Double -> Double
i0' Double
x Double
d Double
ds | Double
ds forall a. Ord a => a -> a -> Bool
< Double
1.0e-30 = Double
1
           | Bool
otherwise = Double
ds forall a. Num a => a -> a -> a
* Double
xforall a. Num a => a -> Int -> a
^!Int
2 forall a. Fractional a => a -> a -> a
/ Double
dforall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
+ (Double -> Double -> Double -> Double
i0' Double
x (Double
dforall a. Num a => a -> a -> a
+Double
2) (Double
ds forall a. Num a => a -> a -> a
* Double
xforall a. Num a => a -> Int -> a
^!Int
2 forall a. Fractional a => a -> a -> a
/ Double
dforall a. Num a => a -> Int -> a
^!Int
2))

-- I don't think this one is correct.  Kay's book uses different variable
-- conventions and I haven't deciphered them yet...

-- | rectangular window

parzen :: Int -- ^ M
       -> Array Int Double -- ^ w[n]

parzen :: Int -> Array Int Double
parzen = (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
parzen'

parzen' :: Double -> Double -> Double
parzen' :: Double -> Double -> Double
parzen' Double
m Double
n =
   if Double
n forall a. Ord a => a -> a -> Bool
<= Double
m forall a. Fractional a => a -> a -> a
/ Double
2
     then Double
2 forall a. Num a => a -> a -> a
* (Double
1forall a. Num a => a -> a -> a
-Double
nforall a. Fractional a => a -> a -> a
/Double
m) forall a. Num a => a -> Int -> a
^! Int
3 forall a. Num a => a -> a -> a
- (Double
1forall a. Num a => a -> a -> a
-Double
2forall a. Num a => a -> a -> a
*Double
nforall a. Fractional a => a -> a -> a
/Double
m) forall a. Num a => a -> Int -> a
^! Int
3
     else Double
2 forall a. Num a => a -> a -> a
* (Double
1forall a. Num a => a -> a -> a
-Double
nforall a. Fractional a => a -> a -> a
/Double
m) forall a. Num a => a -> Int -> a
^! Int
3


makeArray :: (Double -> Double -> Double) -> Int -> Array Int Double
makeArray :: (Double -> Double -> Double) -> Int -> Array Int Double
makeArray Double -> Double -> Double
win Int
m =
   let md :: Double
md = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m
   in  forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
m) forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
win Double
md forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral) [(Int
0::Int) ..]