## Discrete-time LTI systems in Haskell II

Yesterday I wrote a post on how to implement discrete-time LTI systems in Haskell. The code I presented did not quite please my delicate aesthetic sensitivity, as it bundled the declaration of the signals with the declaration of the LTI system under study. This is not practical. Ideally, I would like to define the LTI system in a Haskell script, run it on the GHCi interpreter, then create the signals using the GHCi command line. Why? Because I can then compute the response of the system to various test input signals without reloading the script.

In this post, we will again consider the first-order causal discrete-time LTI system described by the difference equation

$y (n) - \alpha y (n-1) = u (n)$

where $|\alpha| < 1$. In my previous post, I obtained a state-space representation of this LTI system. In this post, I will propose a better implementation of the given system in Haskell.

__________

Signal and System type synonyms

What is a discrete-time signal? We can think of a discrete-time signal as a sequence of numbers or symbols. In Haskell, we often use lists of floating-point numbers to represent real-valued discrete-time signals. More generally, we have the type synonym:

type Signal a = [a]

which says that a signal of type $a$ is a list of elements of type $a$. We can then have signals of various types: integer, fixed-point, floating-point, bit, binary word, etc. This should be particularly useful in case we want to simulate mixed-signal circuits.

What is a discrete-time system? It is that which maps input discrete-time signals to output discrete-time signals. We will use the type synonym:

type System a b = (Signal a -> Signal b)

which tells us that a system of type $a \,\, b$ is a function that maps an input signal of type $a$ to an output signal of type $b$. In other words, it takes a list whose elements are of type $a$, and it returns a list whose elements are of type $b$.

__________

Implementing the LTI system in Haskell

In my previous post, I created the state and output sequences using

xs = scanl f x0 us
ys = zipWith g xs us

Combining the two lines above into a single line, we obtain

ys = zipWith g (scanl f x0 us) us

which returns the output sequence when given the input sequence. Note that the state sequence is not created. Assuming zero initial condition, we can then create a system (i.e., a function that takes a list and returns a list) as follows

sys :: Floating a => System a a
sys us = zipWith g (scanl f 0.0 us) us

which takes signals of type $a$ and returns signals of the same type, where $a$ is a floating-point type (either Float or Double).

Putting it all together in a script, we finally obtain:

type Signal a = [a]

type System a b = (Signal a -> Signal b)

-- state-space model
f,g :: Floating a => a -> a -> a
f x u = (0.5 * x) + u  	-- state-transition
g x u = f x u		-- output

-- define LTI system
sys :: Floating a => System a a
sys us = zipWith g (scanl f 0.0 us) us

which defines the LTI system under study with $\alpha = 1/2$. For this choice of $\alpha$, the system is a lowpass filter.

__________

Example: lowpass-filtering a square wave

We run the script above on GHCi to create the desired LTI system. Let us now create an input signal, a square wave of period equal to 16 and duty cycle equal to 50%, and lowpass-filter it using the LTI system we created:

*Main> -- create square wave
*Main> let ones  = take 8 $repeat 1.0 :: Signal Float *Main> let zeros = take 8$ repeat 0.0 :: Signal Float
*Main> let us = cycle (ones ++ zeros) :: Signal Float
*Main> -- create output sequence
*Main> let ys = sys us
*Main> -- check types
*Main> :type us
us :: Signal Float
*Main> :type ys
ys :: Signal Float
*Main> :type sys
sys :: Floating a => System a a

Finally, we compute the first 64 samples of the input and output signals using function take, and plot the signals using MATLAB:

In case you are wondering why the peak amplitude of the output signal is twice that of the input signal, compute the transfer function of the LTI system and note that the DC gain is equal to $\frac{1}{1 - \alpha} = 2$.