We recently studied the first-order causal LTI system which is described by the difference equation , where . The system can be represented by the block diagram
Observing the block diagram, we conclude that there are three basic operations: addition, multiplication by a constant coefficient, and delay. We arrive at the same conclusion if we rewrite the difference equation in the form . Do note that the output is obtained by adding the input to a scaled and delayed version of the output and, therefore, we have a feedback system.
The system’s input-output relationship can be written as follows
where is the unit-delay operator. Note that we have signals rather than signal samples on both sides of the equation. To clarify, when I say “sample”, I mean the value of the signal at some (discrete) time instant. Let us introduce the linear operator such that the output signal can be written as a function of the input, , assuming zero initial condition. Hence, we obtain . It would be convenient to introduce also a gain operator to carry out the multiplication by a constant coefficient, i.e., . Composing the operators, we finally obtain the following equation
which we will now implement in Haskell.
Implementation in Haskell
Our first implementation of the LTI system under study relied on state-space models and the scanl trick to propagate the initial state forwards in time. Our second implementation was little more than a beautified version of the first one. This third implementation will be radically different from the previous two.
Let us start with the following type synonyms:
type Signal a = [a] type System a b = Signal a -> Signal b
which hopefully will make the code more readable. We build a function that takes two discrete-time signals (of the same type) and returns their elementwise addition:
(.+) :: Num a => Signal a -> Signal a -> Signal a (.+) = zipWith (+)
Since we represent discrete-time signals as lists, the function above merely adds two lists elementwise. Let us test this function:
*Main> -- create test input signals *Main> let us = repeat 1.0 :: Signal Float *Main> let vs = [1..] :: Signal Float *Main> -- add two signals elementwise *Main> let ys = us .+ vs *Main> take 10 ys [2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0] *Main> -- check types *Main> :type (us,vs,ys) (us,vs,ys) :: (Signal Float, Signal Float, Signal Float)
We now implement the unit-delay with zero initial condition:
delay :: Num a => System a a delay us = 0 : us
which right-shifts the input list and introduces a zero at the output list’s head. If we want a unit-delay operator with non-zero initial condition, we would have the following code instead:
delay :: Num a => System a a delay us = ini : us
where ini is the initial condition of the delay block. Lastly, we create the gain operator:
gain :: Num a => a -> System a a gain alpha = map (alpha*)
which takes a number (the gain factor) and returns a system (that maps signals to signals). Note that we use partial function application. One can think of the gain operator as a function that takes a number and a signal and returns a signal. If we fix the first argument (the gain factor), we obtain a function that maps signals to signals, i.e., a system. Here is a quick test of the delay and gain operators:
*Main> -- create signal *Main> let xs = [1..] :: Signal Float *Main> -- delay signal *Main> let ys = delay xs *Main> take 10 ys [0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0] *Main> -- amplify delayed signal *Main> let zs = gain 2.0 ys *Main> take 10 zs [0.0,2.0,4.0,6.0,8.0,10.0,12.0,14.0,16.0,18.0]
Finally, we have the following Haskell script:
type Signal a = [a] type System a b = Signal a -> Signal b -- signal adder (.+) :: Num a => Signal a -> Signal a -> Signal a (.+) = zipWith (+) -- delay operator delay :: Num a => System a a delay us = 0 : us -- gain operator gain :: Num a => a -> System a a gain alpha = map (alpha*) -- build feedback system sys :: Floating a => System a a sys us = us .+ (gain 0.5 . delay . sys) us
where we used . Note that the last line is a direct translation of the equation to Haskell! Beautiful!
To finalize, let us obtain the impulse response of the LTI system under study for :
*Main> -- create unit impulse *Main> let delta = 1.0 : repeat 0.0 :: Signal Float *Main> -- compute impulse response *Main> let hs = sys delta *Main> take 8 hs [1.0,0.5,0.25,0.125,6.25e-2,3.125e-2,1.5625e-2,7.8125e-3]
which is the expected impulse response for . Frankly, I am in awe. It is amazing that this implementation works!