Graph.hs

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

import qualified Prelude as P
import PreludeBase
import NumericPrelude

import Filter.Basic(Filter,apply,transferFunction)
import qualified Data.Map as Map
import Data.Map(Map)
import MathObj.DiscreteMap()  {- Module.C instances for Map -}

import qualified Number.Complex         as Complex
import qualified Algebra.RealField      as RealField
import qualified Algebra.Additive       as Additive
import qualified Algebra.Transcendental as Trans
import qualified Algebra.Module         as Module
import Algebra.Module((*>))
import Orthogonals(Scalar,inverse,add_to_diagonal)


{-|
A filter network is a graph consisting
of nodes (input and output signals)
and edges (filter processes).
Output signals can be taken from every node,
inputs can be injected in some nodes
which means that the signal at this node is superposed with
the injected signal.
The same can be achieved by duplicating the network,
one duplicate per input,
and superposing the corresponding outputs.
It is also sensible to setup a graph without inputs,
e.g. a recursive filter with some initial conditions
that works independent from any input.

In opposition to electric filter networks
digital filter networks must be directed.

Test-case: leap-frog filters like

>     +-----------[d]-----------+
>     v                         |
>(u) -+-> [a] (v) -+-> [b] (w) -+-> [c] (y) -+->
>                  ^                         |
>                  +-----------[e]-----------+

@
v = a·(u + d·w)
w = b·(v + e·y)
y = c· w
@

We model the general network by a list of nodes,
where each node is an adder that holds a list of its inputs.
Each input of a node is an output
of another node that has gone through a processor.
Additionally there may be one input from outside.
In principle a processor could be a simple filter network
as defined by the structure 'Filter'.

The network is an applyable filter
whenever each circle contains a delay.
To compute the transfer function
we have to solve a system of linear equations
which we can construct quite straight forward
from the processors' input lists.

The current design can be abstractly seen
as the system of linear equations:

  y = A*y + u

where A is a matrix containing the edges hosting the filters,
y the vector of output signals,
u the vector of input signals.
In this formulation the number of inputs and outputs must match
but since you are free to ignore some of the inputs and outputs
you can use nodes for pure output, pure input or both of them.

-}

newtype T filter i t a v =
   C (Map i
      [(i,          {- index of the processor whose output goes in here -}
        filter t a v  {- description of the filter -}
        )])


newtype Signal list i v = Signal (Map i (list v))


fromList :: (Ord i) => [(i, [(i, filter t a v)])] -> T filter i t a v
fromList = C . Map.fromList

toList :: T filter i t a v -> [(i, [(i, filter t a v)])]
toList (C fg) = Map.toList fg

signalFromList :: (Ord i) => [(i, list v)] -> Signal list i v
signalFromList = Signal . Map.fromList

signalToList :: Signal list i v -> [(i, list v)]
signalToList (Signal x) = Map.toList x

lookupSignal :: (Ord i) => Signal list i v -> i -> Maybe (list v)
lookupSignal (Signal x) = flip Map.lookup x


{-
  These instance may help to include FilterGraphs
  in even bigger structures.
-}
instance (Ord i, Additive.C (list v), Eq (list v))
            => Additive.C (Signal list i v) where
   zero                      = Signal Additive.zero
   (+) (Signal x) (Signal y) = Signal ((Additive.+) x y)
   (-) (Signal x) (Signal y) = Signal ((Additive.-) x y)
   negate (Signal x)         = Signal (Additive.negate x)

instance (Ord i, Eq a, Additive.C a, Additive.C (list v), Eq (list v),
              Module.C a v, Module.C a (list v))
                 => Module.C a (Signal list i v) where
   s *> (Signal x) = Signal (s *> x)


{-
   It would be interesting to make FilterGraphs
   an instance of Filter.
   To achieve that we had to make GraphSignals an instance of Module.C
   and the transferFunction would no longer return a factor
   but a function that maps input amplitudes
   of a given frequency to output amplitudes.

instance (Ord i, Show i, Filter list filter) =>
      Filter (Signal list i) (T filter i) where
-}

apply :: (Ord i, Show i, Additive.C t, Trans.C t, RealField.C t,
          Module.C a v, Module.C a (list v),
          Filter list filter) =>
   T filter i t a v -> Signal list i v -> Signal list i v
apply (C fg) (Signal inputs) =
   let getInput  i = Map.findWithDefault Additive.zero i inputs
       getOutput i = Map.findWithDefault
                         (error ("Unknown processor: "++show i)) i outputs
       output i edges =
          foldl (Additive.+) (getInput i) (map (\(j,f) ->
                   Filter.Basic.apply f (getOutput j)) edges)
       outputs = Map.mapWithKey output fg
   in  Signal outputs

{-| Compute a matrix that tells how an input frequency
    is mapped to the various output nodes.

   According to the formulation given above
   we have to invert the matrix (I-A).

   Currently this is done by a QR decomposition for each frequency.
   It would be probably faster if we decompose
   the matrix containing polynomial elements.
   Then the inverted matrix would consist of some
   polynomial ratios which can be evaluated for each frequency.
-}
transferFunction ::
   (Ord i, Show i, Trans.C t,
    P.Fractional (Complex.T t), Scalar (Complex.T t),
      Module.C a t, Filter list filter) =>
         T filter i t a v -> t -> [[Complex.T t]]
transferFunction (C fg) w =
   let keys = Map.keys fg
       elts = Map.elems fg
       inputsToMap procs =
          Map.mapWithKey (\_ f -> Filter.Basic.transferFunction f w)
                         (Map.fromList procs)
       makeRow procs =
          map (flip (Map.findWithDefault 0) (inputsToMap procs)) keys
       matrix = map makeRow elts
   in  inverse (add_to_diagonal (-1) matrix)