Composition.hs

{-# OPTIONS -fglasgow-exts -fno-implicit-prelude -fallow-undecidable-instances #-}
module Filter.Composition where

import Filter.Basic (Filter,apply,transferFunction)

import qualified Algebra.Module         as Module
import qualified Algebra.Transcendental as Trans
import qualified Algebra.RealField      as RealField
import qualified Algebra.Field          as Field
import qualified Algebra.Additive       as Additive
import qualified Number.Complex         as Complex

import Algebra.Additive ((+))

import PreludeBase
import NumericPrelude

{- todo:
    - functions that build a FilterComposition for specific filters
        (1st order, universal, allpass, butterworth, chebyshev)
    - functions that turn physical filter parameters into
        internal ones
    - How can these function be combined?
      A function like
         [ FilterComposition v [m] ] -> FilterComposition v [[m]]
      is not satisfying, since the conversion function cannot rely
      that the structure of all FilterComposition v [m] is equal.
      If the list is empty the structure can't even be reconstructed.
-}

{-|
  This describes a generic filter with one input and one main output
  that consists of non-recursive and recursive parts.
  If you use Feedback, make sure that at least
  one of the filters of a circle includes a delay,
  otherwise the recursion will fail.
  The main output is used to glue different parts together.
  Additionally the functions 'apply' and 'transferFunction'
  provide the signals at every node of the network.
-}
data T filter t a v =
     Prim (filter t a v)
       {-^ a filter primitve -}
   | Serial   [T filter t a v]
       {-^ serial chain of filters -}
   | Parallel [T filter t a v]
       {-^ filters working parallel, there output is mixed together -}
   | Feedback (T filter t a v) (T filter t a v)
       {-^ filter the signal in the forward direction and
           feed back the output signal filtered by the second filter -}

{-|
  This is the data structure is used for the results
  of 'apply' and 'transferFunction'.
  Each constructor corresponds to one of 'Filter.Composition.T'.
  By choosing only some of the outputs
  the lazy evaluation will content
  with applying the necessary filter steps, only.
-}
data Sockets s = Sockets {output :: s, socket :: SocketSpec s}

data SocketSpec s =
     Output
   | Multiplier [Sockets s]
   | Adder [Sockets s]
   | Loop (Sockets s) (Sockets s)

instance (Filter list filter) =>
      Filter (list) (T filter) where
{-
   apply :: (Module.C a v) =>
      FilterComposition a v -> TwoWayList v -> TwoWayList v
-}
   apply f x = output (applyMulti f x)
{-
   transferFunction :: (Trans.C b, Module.C a (Complex.T b)) =>
      T filter a v -> b -> (Complex.T b)
-}
   transferFunction f w = output (transferFunctionMulti f w)


{-| Apply a filter network to a signal and keep the output of all nodes.
    Generic function that is wrapped by 'apply'. -}
applyMulti :: (RealField.C t, Trans.C t,
      Module.C a v, Module.C a (list v), Filter list filter) =>
   T filter t a v -> list v -> Sockets (list v)
applyMulti (Prim f) x =
   Sockets (apply f x) Output
applyMulti (Serial fs) x =
   let sq = scanl (\(Sockets y _) -> flip applyMulti y) (Sockets x Output) fs
   in  Sockets (output (last sq)) (Multiplier (tail sq))
applyMulti (Parallel fs) x =
   let socks = map (flip applyMulti x) fs
       y = foldr (Additive.+) zero (map output socks)
   in  Sockets y (Adder socks)
{- the distinction between 'feed' and 'back'
   can be dropped in a more general net structure -}
applyMulti (Feedback feed back) x =
   let sockY@(Sockets y _) = applyMulti feed ((Additive.+) x z)
       sockZ@(Sockets z _) = applyMulti back y
   in  Sockets y (Loop sockY sockZ)


transferFunctionMulti ::
   (Trans.C t, Module.C a t, Filter list filter) =>
      T filter t a v -> t -> Sockets (Complex.T t)
transferFunctionMulti f w = tfAbsolutize 1 (tfRelative w f)

{-| Compute the transitivity for each part of the filter network.
    We must do this in such a relative manner to be able
    to compute feedback. -}
tfRelative ::
   (Trans.C t, Module.C a t, Filter list filter) =>
      t -> T filter t a v -> Sockets (Complex.T t)
tfRelative w (Prim f) =
   Sockets (Filter.Basic.transferFunction f w) Output
tfRelative w (Serial fs) =
   let sq = map (tfRelative w) fs
   in  Sockets (product (map output sq)) (Multiplier sq)
tfRelative w (Parallel fs) =
   let sq = map (tfRelative w) fs
   in  Sockets (sum (map output sq)) (Adder sq)
tfRelative w (Feedback feed back) =
   let sockY = tfRelative w feed
       sockZ = tfRelative w back
       q = output sockY / (1 - output sockZ)
   in  Sockets q (Loop sockY sockZ)


{-| Make the results from 'tfRelative' absolute. -}
tfAbsolutize :: (Field.C a) => a -> Sockets a -> Sockets a
tfAbsolutize x (Sockets y spec) = Sockets (x*y)
   (case spec of
      (Multiplier socks) ->
         let sq = scanl (\(Sockets z _) -> tfAbsolutize z)
                        (Sockets x Output) socks
         in  Multiplier (tail sq)
      (Adder socks) ->
         let sq = map (tfAbsolutize x) socks
         in  Adder sq
      (Loop feed back) ->
         let sockY = tfAbsolutize (x / (1 - output back)) feed
             sockZ = tfAbsolutize (output sockY) back
             -- it should be  x*y == output sockY
         in  Loop sockY sockZ
      Output -> spec)