{-# 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.Allpass where

import qualified Synthesizer.Plain.Signal   as Sig
import qualified Synthesizer.Plain.Modifier as Modifier
import qualified Synthesizer.Causal.Process as Causal
import qualified Synthesizer.Interpolation.Class as Interpol

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 qualified Algebra.Ring                  as Ring
import qualified Algebra.Additive              as Additive

import Algebra.Module((*>))

import Number.Complex ((+:))
import qualified Number.Complex as Complex
import Data.Tuple.HT (mapSnd, )
import Data.Function.HT (nest, )
import Data.List.HT (mapAdjacent, switchR, )

import Control.Monad.Trans.State (State, state, runState, evalState, )

import qualified Prelude as P
import PreludeBase
import NumericPrelude



newtype Parameter a =
   Parameter {getParameter :: a}  {- ^ Feedback factor. -}
   deriving Show


instance Interpol.C a v => Interpol.C a (Parameter v) where
   {-# INLINE scaleAndAccumulate #-}
   scaleAndAccumulate = Interpol.makeMac Parameter getParameter


{-
Shall the phase parameter be of type Phase?
I think no, because for the allpass cascade we divide by the order
and then there is a difference between phase pi and 3*pi.
-}
{-# INLINE parameter #-}
parameter :: Trans.C a =>
     a    {- ^ The phase shift to be achieved for the given frequency in radians. -}
  -> a    {- ^ The frequency we specified the phase shift for. -}
  -> Parameter a
parameter phase frequency =
   let omega = frequency * 2 * pi
       k = (cos phase - cos omega) / (1 - cos (phase - omega))
   in  Parameter k

{-
cos phi = (1-r^2)/(1+r^2)
cos omega = (1-s^2)/(1+s^2)
cos (phi-omega)
   = cos phi * cos omega + sin phi * sin omega
   = ((1-r^2)*(1-s^2) + 4*r*s) / ((1+r^2) * (1+s^2))
k = ((1-r^2)*(1+s^2) - (1+r^2)*(1-s^2)) /
    ((1+r^2) * (1+s^2) - ((1-r^2)*(1-s^2) + 4*r*s))
k = 2*(s^2-r^2) / (2*r^2+2*s^2 - 4*r*s)
k = (s^2-r^2) / (r-s)^2
k = (s+r) / (s-r)
-}
{-# INLINE parameterAlt #-}
parameterAlt :: Trans.C a =>
     a    {- ^ The phase shift to be achieved for the given frequency. -}
  -> a    {- ^ The frequency we specified the phase shift for. -}
  -> Parameter a
parameterAlt phase frequency =
   let s = tan (frequency * pi)
       r = tan (phase/2)
   in  Parameter $ (s+r) / (s-r)

{- |
An approximation to 'parameter' for small phase and frequency values.
It needs only field operations,
however it also needs 'pi', thus the transcendental constraint.
-}
{-# INLINE parameterApprox #-}
parameterApprox :: Trans.C a =>
     a    {- ^ The phase shift to be achieved for the given frequency. -}
  -> a    {- ^ The frequency we specified the phase shift for. -}
  -> Parameter a
parameterApprox phase frequency =
   let omega = frequency * 2 * pi
       k = (omega + phase) / (omega - phase)
   in  Parameter k


-- * atomic first order allpass

{-# INLINE firstOrderStep #-}
firstOrderStep :: (Ring.C a, Module.C a v) =>
   Parameter a -> v -> State (v,v) v
firstOrderStep (Parameter k) u0 =
   state (\(u1,y1) -> let y0 = u1 + k *> (u0-y1) in (y0,(u0,y0)))

{-# INLINE firstOrderModifier #-}
firstOrderModifier :: (Ring.C a, Module.C a v) =>
   Modifier.Simple (v,v) (Parameter a) v v
firstOrderModifier =
   Modifier.Simple (zero,zero) firstOrderStep

{-# INLINE firstOrderCausal #-}
firstOrderCausal :: (Ring.C a, Module.C a v) =>
   Causal.T (Parameter a, v) v
firstOrderCausal =
   Causal.fromSimpleModifier firstOrderModifier

{-# INLINE firstOrder #-}
firstOrder :: (Ring.C a, Module.C a v) =>
   Sig.T (Parameter a) -> Sig.T v -> Sig.T v
firstOrder = Sig.modifyModulated firstOrderModifier


{-# INLINE makePhase #-}
makePhase :: RealTrans.C a => Parameter a -> a -> a
makePhase (Parameter k) frequency =
   let omega = 2*pi * frequency
   in  2 * Complex.phase ((k+cos omega)+:(- sin omega)) + omega


-- * allpass cascade with uniform control

{-# INLINE cascadeParameter #-}
cascadeParameter :: Trans.C a =>
     Int  {- ^ The number of equally designed 1st order allpasses. -}
  -> a    {- ^ The phase shift to be achieved for the given frequency in radians. -}
  -> a    {- ^ The frequency we specified the phase shift for. -}
  -> Parameter a
cascadeParameter order phase =
   parameter (phase / fromIntegral order)

{-# INLINE cascadeParameterAlt #-}
cascadeParameterAlt :: Trans.C a =>
     Int  {- ^ The number of equally designed 1st order allpasses. -}
  -> a    {- ^ The phase shift to be achieved for the given frequency in radians. -}
  -> a    {- ^ The frequency we specified the phase shift for. -}
  -> Parameter a
cascadeParameterAlt order phase frequency =
   let orderFloat = fromIntegral order
       omega = frequency * 2 * pi
       phi = phase / orderFloat
       k = (cos phi - cos omega) / (1 - cos (phi - omega))
   in  Parameter k

{-# INLINE flangerPhase #-}
flangerPhase :: Trans.C a => a
flangerPhase = -2*pi

{-# INLINE flangerParameter #-}
flangerParameter :: Trans.C a => Int -> a -> Parameter a
flangerParameter order frequency =
   cascadeParameter order flangerPhase frequency


{-# INLINE cascadeStep #-}
cascadeStep :: (Ring.C a, Module.C a v) =>
   Parameter a -> v -> State [v] v
cascadeStep = cascadeStepRec

{-
internal storage is not very efficient
because the second value of one pair is equal
to the first value of the subsequent value
-}
{-# INLINE cascadeStepStackPairs #-}
cascadeStepStackPairs :: (Ring.C a, Module.C a v) =>
   Parameter a -> v -> State [(v,v)] v
cascadeStepStackPairs k =
   -- stackStatesR would work as well, but with reversed list of states
   Modifier.stackStatesL (firstOrderStep k)

{-# INLINE cascadeStepStack #-}
{-# INLINE cascadeStepRec #-}
{-# INLINE cascadeStepScanl #-}
cascadeStepStack, cascadeStepRec, cascadeStepScanl ::
   (Ring.C a, Module.C a v) =>
   Parameter a -> v -> State [v] v
cascadeStepStack k x =
   state $
      mapSnd fromPairs .
      runState (cascadeStepStackPairs k x) .
      toPairs

{-# INLINE fromPairs #-}
fromPairs :: [(a,a)] -> [a]
fromPairs xs@(x:_) = fst x : map snd xs
fromPairs [] = error "Allpass.fromPairs: empty list"

{-# INLINE toPairs #-}
toPairs :: [a] -> [(a,a)]
toPairs xs = mapAdjacent (,) xs


cascadeStepRec (Parameter k) x = state $ \s ->
   let crawl _ [] = error "Allpass.crawl needs at least one element in the list"
       crawl u0 (_:[]) = u0:[]
       crawl u0 (u1:y1:us) =
           let y0 = u1 + k *> (u0-y1)
           in  u0 : crawl y0 (y1:us)
       news = crawl x s
   in  (last news, news)

cascadeStepScanl k x = state $ \s ->
   let news =
          scanl
             (evalState . firstOrderStep k)
             x (mapAdjacent (,) s)
   in  (switchR
           (error "Allpass.cascade needs at least one element in the state list")
           (flip const) news,
        news)


{-# INLINE cascadeModifier #-}
cascadeModifier :: (Ring.C a, Module.C a v) =>
   Int -> Modifier.Simple [v] (Parameter a) v v
cascadeModifier order =
   Modifier.Simple (replicate (succ order) zero) cascadeStep

{-# INLINE cascadeCausal #-}
{-# INLINE cascadeCausalStacked #-}
{-# INLINE cascadeCausalModifier #-}
cascadeCausal, cascadeCausalStacked, cascadeCausalModifier ::
   (Ring.C a, Module.C a v) =>
   Int -> Causal.T (Parameter a, v) v
cascadeCausal = cascadeCausalModifier

cascadeCausalStacked order =
   Causal.replicateControlled order firstOrderCausal

cascadeCausalModifier order =
   Causal.fromSimpleModifier (cascadeModifier order)


{-# INLINE cascade #-}
{-# INLINE cascadeState #-}
{-# INLINE cascadeIterative #-}
cascade, cascadeState, cascadeIterative ::
   (Ring.C a, Module.C a v) =>
   Int -> Sig.T (Parameter a) -> Sig.T v -> Sig.T v

{-| Choose one of the implementations below -}
cascade = cascadeState

{-| Simulate the Allpass cascade by a list of states of the partial allpasses -}
cascadeState order =
   Sig.modifyModulated (cascadeModifier order)

{-| Directly implement the allpass cascade as multiple application of allpasses of 1st order -}
cascadeIterative order c =
   nest order (firstOrder c)



-- * allpass cascade with independently controlled atomic allpasses

{-# INLINE cascadeDiverseStep #-}
{-# INLINE cascadeDiverseStepScanl #-}
cascadeDiverseStep, cascadeDiverseStepScanl :: (Ring.C a, Module.C a v) =>
   [Parameter a] -> v -> State [v] v
cascadeDiverseStep = cascadeDiverseStepScanl

cascadeDiverseStepScanl ks x = state $ \s ->
   let news =
          scanl
             (\u0 (k,uy1) -> evalState (firstOrderStep k u0) uy1)
             x (zip ks $ mapAdjacent (,) s)
   in  (switchR
           (error "Allpass.cascadeDiverse needs at least one element in the state list")
           (flip const) news,
        news)