-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.Goertzel
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This is an implementation of Goertzel's algorithm, which computes one
-- bin of a DFT.  A description can be found in Oppenheim and Schafer's
-- /Discrete Time Signal Processing/, pp 585-587.
--
-----------------------------------------------------------------------------

-- TODO: do the cipherin' to figure out the best simplification for the
-- cgoertzel_power case

-- TODO: Bonzanigo's phase correction

module Numeric.Transform.Fourier.Goertzel where

import Data.Array
import Data.Complex

-- | Goertzel's algorithm for complex inputs

cgoertzel :: (RealFloat a, Ix b, Integral b) => Array b (Complex a) -- ^ x[n]
	  -> b -- ^ k
	  -> Complex a -- ^ X[k]

cgoertzel :: forall a b.
(RealFloat a, Ix b, Integral b) =>
Array b (Complex a) -> b -> Complex a
cgoertzel Array b (Complex a)
x0 b
k = [Complex a] -> Complex a -> Complex a -> Complex a
g (forall i e. Array i e -> [e]
elems Array b (Complex a)
x0) Complex a
0 Complex a
0
    where w :: a
w = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral b
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
          a :: a
a = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos a
w
	  g :: [Complex a] -> Complex a -> Complex a -> Complex a
g []     Complex a
x1 Complex a
x2 = Complex a
x1 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> Complex a
cis a
w forall a. Num a => a -> a -> a
- Complex a
x2
	  g (Complex a
x:[Complex a]
xs) x1 :: Complex a
x1@(a
x1r:+a
x1i) Complex a
x2 = [Complex a] -> Complex a -> Complex a -> Complex a
g [Complex a]
xs (Complex a
x forall a. Num a => a -> a -> a
+ (a
aforall a. Num a => a -> a -> a
*a
x1rforall a. a -> a -> Complex a
:+a
aforall a. Num a => a -> a -> a
*a
x1i) forall a. Num a => a -> a -> a
- Complex a
x2) Complex a
x1
	  n :: b
n = (forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array b (Complex a)
x0) forall a. Num a => a -> a -> a
- b
1

-- | Power via Goertzel's algorithm for complex inputs

cgoertzel_power :: (RealFloat a, Ix b, Integral b) => Array b (Complex a) -- ^ x[n]
		-> b -- ^ k
		-> a -- ^ |X[k]|^2

cgoertzel_power :: forall a b.
(RealFloat a, Ix b, Integral b) =>
Array b (Complex a) -> b -> a
cgoertzel_power Array b (Complex a)
x b
k = (forall a. RealFloat a => Complex a -> a
magnitude forall a b. (a -> b) -> a -> b
$ forall a b.
(RealFloat a, Ix b, Integral b) =>
Array b (Complex a) -> b -> Complex a
cgoertzel Array b (Complex a)
x b
k)forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int)

-- | Goertzel's algorithm for real inputs

rgoertzel :: (RealFloat a, Ix b, Integral b) => Array b a -- ^ x[n]
	  -> b -- ^ k
	  -> Complex a -- ^ X[k]

rgoertzel :: forall a b.
(RealFloat a, Ix b, Integral b) =>
Array b a -> b -> Complex a
rgoertzel Array b a
x0 b
k = [a] -> a -> a -> Complex a
g (forall i e. Array i e -> [e]
elems Array b a
x0) a
0 a
0
    where w :: a
w = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral b
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
          a :: a
a = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos a
w
	  g :: [a] -> a -> a -> Complex a
g []     a
x1 a
x2 = ((a
x1 forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
cos a
w forall a. Num a => a -> a -> a
* a
x2) forall a. a -> a -> Complex a
:+ a
x2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin a
w)
	  g (a
x:[a]
xs) a
x1 a
x2 = [a] -> a -> a -> Complex a
g [a]
xs (a
x forall a. Num a => a -> a -> a
+ a
a forall a. Num a => a -> a -> a
* a
x1 forall a. Num a => a -> a -> a
- a
x2) a
x1
	  n :: b
n = (forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array b a
x0) forall a. Num a => a -> a -> a
- b
1

-- | Power via Goertzel's algorithm for real inputs

rgoertzel_power :: (RealFloat a, Ix b, Integral b) => Array b a -- ^ x[n]
		-> b -- ^ k
		-> a -- ^ |X[k]|^2

rgoertzel_power :: forall a b. (RealFloat a, Ix b, Integral b) => Array b a -> b -> a
rgoertzel_power Array b a
x0 b
k = [a] -> a -> a -> a
g (forall i e. Array i e -> [e]
elems Array b a
x0) a
0 a
0
    where w :: a
w = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral b
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
          a :: a
a = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos a
w
	  g :: [a] -> a -> a -> a
g []     a
x1 a
x2 = a
x1forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) forall a. Num a => a -> a -> a
+ a
x2forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) forall a. Num a => a -> a -> a
- a
a forall a. Num a => a -> a -> a
* a
x1 forall a. Num a => a -> a -> a
* a
x2
	  g (a
x:[a]
xs) a
x1 a
x2 = [a] -> a -> a -> a
g [a]
xs (a
x forall a. Num a => a -> a -> a
+ a
a forall a. Num a => a -> a -> a
* a
x1 forall a. Num a => a -> a -> a
- a
x2) a
x1
	  n :: b
n = (forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array b a
x0) forall a. Num a => a -> a -> a
- b
1