We recently studied the first-order causal LTI system which is described by the difference equation $y (n) - \alpha y (n-1) = u (n)$, where $|\alpha| < 1$. 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 $y (n) = u (n) + \alpha y (n-1)$. 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

$y = u + \alpha\,\mathcal{D} (y)$

where $\mathcal{D}$ 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 $\mathcal{H}$ such that the output signal can be written as a function of the input, $y = \mathcal{H} (u)$, assuming zero initial condition. Hence, we obtain $\mathcal{H} (u) = u + \alpha\,\mathcal{D} (\mathcal{H} (u))$. It would be convenient to introduce also a gain operator $\mathcal{G}_{\alpha}$ to carry out the multiplication by a constant coefficient, i.e., $\mathcal{G}_{\alpha} (x) = \alpha \, x$. Composing the operators, we finally obtain the following equation

$\mathcal{H} (u) = u + (\mathcal{G}_{\alpha} \circ \mathcal{D} \circ \mathcal{H}) (u)$

which we will now implement 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.

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

(.+) :: 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 $\alpha = 0.5$. Note that the last line is a direct translation of the equation $\mathcal{H} (u) = u + (\mathcal{G}_{\alpha} \circ \mathcal{D} \circ \mathcal{H}) (u)$ to Haskell! Beautiful!

To finalize, let us obtain the impulse response of the LTI system under study for $\alpha = 0.5$:

*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 $h (n) = \alpha^n$ for $n \geq 0$. Frankly, I am in awe. It is amazing that this implementation works!