-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Source.Oscillator
-- Copyright   :  (c) Matthew Donadio 1998,2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- NCO and NCOM functions
--
-----------------------------------------------------------------------------

module DSP.Source.Oscillator (nco, ncom,
			      quadrature_nco, complex_ncom,
			      quadrature_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 :: forall a. RealFloat a => a -> a -> [a]
nco a
wn a
phi = [a]
y
    where a0 :: a
a0 = a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos a
wn
	  y1 :: [a]
y1 = -(forall a. Floating a => a -> a
sin (a
wn forall a. Num a => a -> a -> a
+ a
phi)) forall a. a -> [a] -> [a]
: [a]
y
          y2 :: [a]
y2 = -(forall a. Floating a => a -> a
sin (a
2 forall a. Num a => a -> a -> a
* a
wn forall a. Num a => a -> a -> a
+ a
phi)) forall a. a -> [a] -> [a]
: [a]
y1
          y :: [a]
y  = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (forall a b. (a -> b) -> [a] -> [b]
map (a
a0 forall a. Num a => a -> a -> a
*) [a]
y1) [a]
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 :: forall a. RealFloat a => a -> a -> [a] -> [a]
ncom a
wn a
phi [a]
x = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [a]
x (forall a. RealFloat a => a -> a -> [a]
nco a
wn a
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 :: forall a. RealFloat a => Complex a -> Complex a
agc (a
x:+a
y) = a
x forall a. Num a => a -> a -> a
* a
r forall a. a -> a -> Complex a
:+ a
y forall a. Num a => a -> a -> a
* a
r
    where r :: a
r = (a
3 forall a. Num a => a -> a -> a
- a
x forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
- a
y forall a. Num a => a -> a -> a
* a
y) forall a. Fractional a => a -> a -> a
/ a
2

-- | 'quadrature_nco' returns an infinite list representing a complex phasor
-- with a phase step of wn radians, ie a quadrature nco with normalized
-- frequency wn radians\/sample.  Since Haskell uses lazy evaluation,
-- 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 :: forall a. RealFloat a => a -> a -> [Complex a]
quadrature_nco a
wn a
phi = (forall a. Floating a => a -> Complex a
cis a
phi) forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
(*) (forall a. Floating a => a -> Complex a
cis a
wn)) (forall a. RealFloat a => a -> a -> [Complex a]
quadrature_nco a
wn a
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 :: forall a. RealFloat a => a -> a -> [Complex a] -> [Complex a]
complex_ncom a
_  a
_   [] = []
complex_ncom a
wn a
phi [Complex a]
x  = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (forall a. RealFloat a => a -> a -> [Complex a]
quadrature_nco a
wn a
phi) [Complex a]
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 :: forall a. RealFloat a => Complex a -> Complex a -> a
quadrature_mult (a
x1:+a
y1) (a
x2:+a
y2) = a
x1 forall a. Num a => a -> a -> a
* a
x2 forall a. Num a => a -> a -> a
+ a
y1 forall a. Num a => a -> a -> a
* a
y2

-- | 'quadrature_ncom' mixes the complex input x with a quadrature nco with
-- normalized frequency wn radians\/sample in quadrature (I\/Q modulation)

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

quadrature_ncom :: forall a. RealFloat a => a -> a -> [Complex a] -> [a]
quadrature_ncom a
wn a
phi [Complex a]
x = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. RealFloat a => Complex a -> Complex a -> a
quadrature_mult [Complex a]
x (forall a. RealFloat a => a -> a -> [Complex a]
quadrature_nco a
wn a
phi)