{-# LANGUAGE NoImplicitPrelude #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE FlexibleContexts #-} module Synthesizer.Filter.Graph where import qualified Prelude as P import PreludeBase import NumericPrelude import qualified Synthesizer.Filter.Basic as FilterBasic import Synthesizer.Filter.Basic (Filter, apply, ) 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) -> FilterBasic.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 -> FilterBasic.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)