Discrete-time LTI systems in Haskell

Consider the first-order causal discrete-time LTI system described by the following difference equation

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

where $|\alpha| < 1$ . This system can be represented by the block diagram

where $D$ is a unit-delay block, the output of which is $y (n-1)$ . We now introduce state variable $x (n)$ denoting the output of the unit-delay block, i.e., $x(n) = y(n-1)$ . Thus, $x (n+1) = y(n)$ and, since $y (n) = \alpha y (n-1) + u (n)$ , we obtain the following state-transition equation

$x (n+1) = \alpha x (n) + u (n)$ .

The initial condition is $x(0) = y(-1)$ , which we assume to be zero. Note that the output is $y (n) = x(n+1)$ and, thus, the output equation is

$y (n) = \alpha x (n) + u (n)$ .

Lastly, we define $f (x,u) = g (x,u) := \alpha x + u$ , which allows us to write the state-transition and output equations as follows

$\begin{array}{rl} x (n+1) &= f (x(n), u(n))\\ y (n) &= g (x(n), u(n))\\\end{array}$

We can now implement this state-space model in Haskell!

__________

Implementation in Haskell

Given an initial state $x_0$ and an input sequence $u = [u_0, u_1, \dots]$ , we can compute the state sequence using the scanl trick I wrote about last weekend. Once we have obtained the state sequence $x = [x_0, x_1, \dots]$ , the output sequence is

$y = [y_0, y_1, \dots] = [g (x_0, u_0), g (x_1, u_1), \dots]$ .

Do you recognize the pattern? We are zipping lists $x$ and $u$ using function $g$ !! We thus use (higher-order) function zipWith to compute the output sequence.

For $\alpha = 7/8$ and a constant input sequence, the following Haskell code computes the state and output trajectories:

-- constant input sequence
us :: Floating a => [a]
us = repeat 1.0

-- parameter alpha
alpha :: Floating a => a
alpha = 7/8

-- state-transition function
f :: Floating a => a -> a -> a
f x u = (alpha * x) + u

-- output function
g :: Floating a => a -> a -> a
g x u = f x u

-- initial condition
x0 :: Floating a => a
x0 = 0.0

-- state sequence
xs :: Floating a => [a]
xs = scanl f x0 us

-- output sequence
ys :: Floating a => [a]
ys = zipWith g xs us

Not quite. I am still stuck in the imperative paradigm… and need to be corrected. The code above does not compute anything. It does declare what the input, state and output sequences are. After running the script above in GHCi, we compute 30 samples of the input and output sequences using function take:

*Main> take 30 us
[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
*Main> take 30 ys
[1.0,1.875,2.640625,3.310546875,3.896728515625,
4.409637451171875,4.858432769775391,5.251128673553467,
5.5947375893592834,5.895395390689373,6.158470966853201,
6.388662095996551,6.590079333996982,6.7663194172473595,
6.92052949009144,7.05546330383001,7.173530390851258,
7.276839091994852,7.367234205495495,7.446329929808559,
7.515538688582489,7.576096352509678,7.629084308445968,
7.675448769890222,7.716017673653944,7.751515464447201,
7.782576031391301,7.809754027467388,7.833534774033964,
7.854342927279719]

We can plot these sequences in MATLAB:

We now know how to simulate discrete-time LTI systems using Haskell. Note that our framework is based on state-space models, not on transfer functions. Hence, it should be possible to simulate discrete-time nonlinear systems in the same manner: using scanl to compute the state sequence, and zipWith to compute the output sequence. All we need is a state-space representation of the system under study.

Like this:

One blogger likes this post.

Tags: Haskell, List Processing, LTI Systems, Signals & Systems, State-Space Models, Systems Theory

This entry was posted on January 21, 2012 at 04:57 and is filed under Haskell, LTI Systems, Systems Theory. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

Rod Carvalho

Leave a Reply Cancel reply

Follow “Rod Carvalho”