```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Transform.Fourier.Goertzel
--
-- 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 x0 k = g (elems x0) 0 0
where w = 2 * pi * fromIntegral k / fromIntegral n
a = 2 * cos w
g []     x1 x2 = x1 * cis w - x2
g (x:xs) x1@(x1r:+x1i) x2 = g xs (x + (a*x1r:+a*x1i) - x2) x1
n = (snd \$ bounds x0) - 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 x k = (magnitude \$ cgoertzel x k)^(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 x0 k = g (elems x0) 0 0
where w = 2 * pi * fromIntegral k / fromIntegral n
a = 2 * cos w
g []     x1 x2 = ((x1 - cos w * x2) :+ x2 * sin w)
g (x:xs) x1 x2 = g xs (x + a * x1 - x2) x1
n = (snd \$ bounds x0) - 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 x0 k = g (elems x0) 0 0
where w = 2 * pi * fromIntegral k / fromIntegral n
a = 2 * cos w
g []     x1 x2 = x1^(2::Int) + x2^(2::Int) - a * x1 * x2
g (x:xs) x1 x2 = g xs (x + a * x1 - x2) x1
n = (snd \$ bounds x0) - 1
```