{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright   :  (c) Henning Thielemann 2008
License     :  GPL

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes
-}
module Synthesizer.Plain.Filter.Recursive.AllpassPoly where

-- import qualified Synthesizer.Plain.Signal   as Sig
-- import qualified Synthesizer.Plain.Modifier as Modifier

import qualified Algebra.Module                as Module
import qualified Algebra.RealTranscendental    as RealTrans
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.Field                 as Field

import Number.Complex (cis,(+:),real,imag)
import qualified Number.Complex as Complex
import Orthogonals(Scalar,one_ket_solution)

import qualified Prelude as P
import NumericPrelude.Base
import NumericPrelude.Numeric



newtype Parameter a = Parameter [a]
   deriving Show

{- | Compute coefficients for an allpass that shifts low frequencies
     by approximately the shift you want.
     To achieve this we solve a linear least squares problem,
     where low frequencies are more weighted than high ones.
     The output is a list of coefficients for an arbitrary order allpass. -}
shiftParam :: (Scalar a, P.Fractional a, Trans.C a) =>
   Int -> a -> a -> Parameter a
shiftParam order weight phase =
    let {- construct matrix for normal equations -}
        normalVector = map negate
           (map (scalarProdScrewExp weight order phase 0) [1..order])
        normalMatrix = map (\j ->
            map (scalarProdScrewExp weight order phase j) [1..order]) [1..order]
    in  Parameter (one_ket_solution normalMatrix normalVector)

{-
  GNUPlot.plotFunc (GNUPlot.linearScale 500 (0,1)) ((fwrap (-pi,pi)).(makePhase (shiftParam 6 (-6) (-pi/2::Double))))
-}
makePhase :: RealTrans.C a => Parameter a -> a -> a
makePhase (Parameter ks) frequency =
    let omega  = 2*pi * frequency
        omegas = iterate (omega+) omega
        denom = 1+sum (zipWith (\k w -> k*cos w +: k*sin w) ks omegas)
    in  2 * Complex.phase denom - omega*(fromIntegral (length ks))

{- integrate (0,2*pi) (\omega -> exp (r*omega) * screwProd order phase k j omega) -}
scalarProdScrewExp :: Trans.C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp r order phase k j =
    let (intCos,intSin) = integrateScrewExp r (k+j-order)
    in  2 * (fst (integrateScrewExp r (k-j)) -
              (cos phase * intCos + sin phase * intSin))

screwProd :: Trans.C a => Int -> a -> Int -> Int -> a -> a
screwProd order phase k j omega =
    let z0 = cis (fromIntegral k * omega) -
                       cis phase * cis (fromIntegral (order-k) * omega)
        z1 = cis (fromIntegral j * omega) -
                       cis phase * cis (fromIntegral (order-j) * omega)
    in  real z0 * real z1 + imag z0 * imag z1

{- integrate (0,2*pi) (\omega -> (exp (r*omega) +: 0) * cis (k*omega)) -}
integrateScrewExp :: Trans.C a => a -> Int -> (a,a)
integrateScrewExp r kInt =
    let k = fromIntegral kInt
        q = (exp (2*pi*r) - 1) / (r^2 + k^2)
    in  (r*q, -k*q)

{- Should be moved to NumericPrelude.Numeric -}
integrateNum :: (Field.C a, Module.C a v) => Int -> (a,a) -> (a->v) -> v
integrateNum n (lo,hi) f =
    let xs = map (\k -> lo + (hi-lo) * fromIntegral k / fromIntegral n)
                 [1..(n-1)]
    in  ((hi-lo) / fromIntegral n) *>
        (foldl (+) ((1/2 `asTypeOf` lo) *> (f lo + f hi))
               (map f xs))