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