```-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Source.Oscillator
--
-- Stability   :  experimental
-- Portability :  portable
--
-- NCO and NCOM functions
--
-----------------------------------------------------------------------------

module DSP.Source.Oscillator (nco, ncom,
agc) where

import Data.Complex

-- | 'nco' creates a sine wave with normalized frequency wn (numerically
-- controlled oscillator, or NCO) using the recurrence relation y[n] =
-- 2cos(wn)*y[n-1] - y[n-2].  Eventually, cumulative errors will creep
-- into the data.  This is unavoidable since performing AGC on this type
-- of real data is hard.  The good news is that the error is small with
-- floating point data.

nco :: RealFloat a => a -- ^ w
-> a -- ^ phi
-> [a] -- ^ y

nco wn phi = y
where a0 = 2 * cos wn
y1 = -(sin (wn + phi)) : y
y2 = -(sin (2 * wn + phi)) : y1
y  = zipWith (-) (map (a0 *) y1) y2

-- | 'ncom' mixes (multiplies) x by a real sine wave with normalized
-- frequency wn.  This is usually called an NCOM: Numerically Controlled
-- Oscillator and Modulator.

ncom :: RealFloat a => a -- ^ w
-> a -- ^ phi
-> [a] -- ^ x
-> [a] -- ^ y

ncom wn phi x = zipWith (*) x (nco wn phi)

-- agc is used in quadrature_nco (below) to scale a complex phasor to
-- have length as close to 1 as possible, ie perform some automatic gain
-- control.  Since we aren't computing sin and cos for each sample, not
-- using AGC would results in cumulative errors (small one with floating
-- point data).  The Complex class includes the signum function which
-- will do what we want, but we will use the approximation 1/sqrt(x) ~=
-- (3-x)/2 for x ~= 1 to eliminate doing a sqrt for every point.

agc         :: RealFloat a => Complex a -> Complex a
agc (x:+y) = x * r :+ y * r
where r = (3 - x * x - y * y) / 2

-- | 'quadrature_nco' returns an infinite list representing a complex phasor
-- with a phase step of wn radians, ie a quadrature nco with normalized
-- rotate will only be computed once, so this NCO uses only one sin and
-- one cos for the entire list, at the expense of 4 mults, 1 add, and 1
-- subtract per point.

quadrature_nco :: RealFloat a => a -- ^ w
-> a -- ^ phi
-> [ Complex a ] -- ^ y

quadrature_nco wn phi = (cis phi) : map ((*) (cis wn)) (quadrature_nco wn phi)

-- | 'complex_ncom' mixes the complex input x with a quardatue nco with
-- normalized frequency wn radians\/sample using complex multiplies
-- (perform a complex spectral shift)

complex_ncom      :: RealFloat a => a -- ^ w
-> a -- ^ phi
-> [ Complex a ] -- ^ x
-> [ Complex a ] -- ^ y

complex_ncom _  _   [] = []
complex_ncom wn phi x  = zipWith (*) (quadrature_nco wn phi) x

-- quadrature_mults returns the sum of the real parts and the imagimary
-- parts of two complex numbers (dot product)

quadrature_mult                 :: RealFloat a => Complex a -> Complex a -> a
quadrature_mult (x1:+y1) (x2:+y2) = x1 * x2 + y1 * y2

-- | 'quadrature_ncom' mixes the complex input x with a quadrature nco with