Posts Tagged ‘List Processing’

Generating all words over an alphabet

October 1, 2012

An alphabet \Sigma is a finite set of symbols. A word w of length n over an alphabet \Sigma is a sequence of symbols from \Sigma, i.e.,

w : \{0,1,\dots,n-1\} \to \Sigma

Alternatively, we can view a word w of length n over the alphabet \Sigma as an n-tuple (i.e., an ordered list with n elements) over \Sigma, i.e., w \in \Sigma^n. The set of all finite words over \Sigma, including the empty word \varepsilon, is denoted by \Sigma^*, where * is the Kleene star.

We thus encounter the following problem:

Problem: Given an alphabet \Sigma, how do we generate all the words over \Sigma? In other words, given \Sigma, how do we generate the Kleene closure \Sigma^*?

The following Haskell script solves this problem:

g :: [a] -> [[a]] -> [[a]]
g alphabet = concat . map (\xs -> [ xs ++ [s] | s <- alphabet])

allwords :: [a] -> [[a]]
allwords alphabet = concat $ iterate (g alphabet) [[]]

where alphabet must be a finite list, otherwise the execution of the script won’t ever terminate. The script above, although short, uses a lot of machinery: list comprehension, higher-order functions map and iterate, concatenation of lists of lists, and partial function application. Let us test this script. First, we load it into GHCi. Then we can run a GHCi session:

*Main> -- define alphabet
*Main> let alphabet = [0,1]
*Main> -- define function f
*Main> let f = g alphabet
*Main> -- check type
*Main> :t f
f :: [[Integer]] -> [[Integer]]
*Main> -- apply f to several lists of lists
*Main> f [[]]
[[0],[1]]
*Main> (f . f) [[]]
[[0,0],[0,1],[1,0],[1,1]]
*Main> (f . f . f) [[]]
[[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]]

So far, so good. Suppose that we would like to find all words over alphabet \Sigma_2 := \{0,1\} whose length is less than or equal to 4. How many such words are there? There are |\Sigma_2|^k = 2^k words of length k. Hence, there are

\displaystyle\sum_{k=0}^n 2^k = \frac{2^{n+1} -1}{2-1} = 2^{n+1} - 1

binary words of length less than or equal to n. Hence, if we make n = 4, then we obtain 2^{n+1}-1 = 31. Continuing our GHCi session:

*Main> take 31 $ allwords alphabet
[[],[0],[1],
[0,0],[0,1],[1,0],[1,1],
[0,0,0],[0,0,1],[0,1,0],[0,1,1],
[1,0,0],[1,0,1],[1,1,0],[1,1,1],
[0,0,0,0],[0,0,0,1],[0,0,1,0],[0,0,1,1],
[0,1,0,0],[0,1,0,1],[0,1,1,0],[0,1,1,1],
[1,0,0,0],[1,0,0,1],[1,0,1,0],[1,0,1,1],
[1,1,0,0],[1,1,0,1],[1,1,1,0],[1,1,1,1]]

Very nice! We could use this in Coding Theory! For example, suppose that we would like to find all binary words of length 8 whose Hamming weight is equal to 5. There are

\displaystyle\binom{8}{5} = \frac{8!}{3! 5!} = 56

such words. We can find them as follows:

*Main> -- take all words of length less than or equal to 8
*Main>  let words = take 511 $ allwords alphabet
*Main> -- filter out binary words of length 8
*Main> let wordsL8 = filter (\xs -> length xs == 8) words
*Main> length wordsL8
256
*Main> -- filter out binary words of length 8
*Main> -- and also of Hamming weight equal to 5
*Main> let wordsL8H5 = filter (\xs -> sum xs == 5) wordsL8
*Main> length wordsL8H5
56
*Main> wordsL8H5
[[0,0,0,1,1,1,1,1],[0,0,1,0,1,1,1,1],[0,0,1,1,0,1,1,1],
 [0,0,1,1,1,0,1,1],[0,0,1,1,1,1,0,1],[0,0,1,1,1,1,1,0],
 [0,1,0,0,1,1,1,1],[0,1,0,1,0,1,1,1],[0,1,0,1,1,0,1,1],
 [0,1,0,1,1,1,0,1],[0,1,0,1,1,1,1,0],[0,1,1,0,0,1,1,1],
 [0,1,1,0,1,0,1,1],[0,1,1,0,1,1,0,1],[0,1,1,0,1,1,1,0],
 [0,1,1,1,0,0,1,1],[0,1,1,1,0,1,0,1],[0,1,1,1,0,1,1,0],
 [0,1,1,1,1,0,0,1],[0,1,1,1,1,0,1,0],[0,1,1,1,1,1,0,0],
 [1,0,0,0,1,1,1,1],[1,0,0,1,0,1,1,1],[1,0,0,1,1,0,1,1],
 [1,0,0,1,1,1,0,1],[1,0,0,1,1,1,1,0],[1,0,1,0,0,1,1,1],
 [1,0,1,0,1,0,1,1],[1,0,1,0,1,1,0,1],[1,0,1,0,1,1,1,0],
 [1,0,1,1,0,0,1,1],[1,0,1,1,0,1,0,1],[1,0,1,1,0,1,1,0],
 [1,0,1,1,1,0,0,1],[1,0,1,1,1,0,1,0],[1,0,1,1,1,1,0,0],
 [1,1,0,0,0,1,1,1],[1,1,0,0,1,0,1,1],[1,1,0,0,1,1,0,1],
 [1,1,0,0,1,1,1,0],[1,1,0,1,0,0,1,1],[1,1,0,1,0,1,0,1],
 [1,1,0,1,0,1,1,0],[1,1,0,1,1,0,0,1],[1,1,0,1,1,0,1,0],
 [1,1,0,1,1,1,0,0],[1,1,1,0,0,0,1,1],[1,1,1,0,0,1,0,1],
 [1,1,1,0,0,1,1,0],[1,1,1,0,1,0,0,1],[1,1,1,0,1,0,1,0],
 [1,1,1,0,1,1,0,0],[1,1,1,1,0,0,0,1],[1,1,1,1,0,0,1,0],
 [1,1,1,1,0,1,0,0],[1,1,1,1,1,0,0,0]]

where we used function take (again) and higher-order function filter. The two predicates were built using anonymous functions.

I am happy with this script. Please let me know in case you are not.

__________

Related:

My implementation of (!!)

September 20, 2012

In Haskell we can easily create a list and then access its elements using the (!!) function, which is defined in the Prelude. Here is a very brief GHCi session:

Prelude> let xs = [7,8,9]
Prelude> xs !! 0
7
Prelude> xs !! 1
8
Prelude> xs !! 2
9

So far, so good. What if the index is negative or equals / exceeds the list’s length? Let’s see what happens in those cases:

Prelude> xs !! (-1)
*** Exception: Prelude.(!!): negative index

Prelude> xs !! 3
*** Exception: Prelude.(!!): index too large

As expected, we get error messages. What if we used the Maybe data type to avoid exceptions? This is exercise 4 in chapter 3 of O’Donnell & Hall & Page [1], which is phrased as follows:

Write (!!), a function that takes a natural number n and a list and selects the nth element of the list. List elements are indexed from 0, not 1, and since the type of the incoming number does not prevent it from being out of range, the result should be a Maybe type.

The aforementioned authors propose the following implementation:

import Prelude hiding ((!!))

(!!) :: Int -> [a] -> Maybe a
(!!) n [] = Nothing
(!!) 0 (x:xs) = Just x
(!!) n (x:xs) = (!!) (n-1) xs

where I added the import line to hide the standard (!!) function that is defined in the Prelude. My first thought was that the authors switched the function arguments, which makes the function look silly. Let’s give it a try. Here’s another GHCi session:

*Main> let xs = [7,8,9]
*Main> 0 !! xs
Just 7
*Main> 1 !! xs
Just 8
*Main> 2 !! xs
Just 9
*Main> (-1) !! xs
Nothing
*Main> 3 !! xs
Nothing

It appears to be working, but specifying the index before the list looks rather ugly. Wait, what if the index is negative? For example, why does (-3) !! xs return Nothing? Let’s use equational reasoning to find out:

(-3) !! [7,8,9] = (-4) !! [8,9] = 
                = (-5) !! [9] =
                = (-6) !! [] = 
                = Nothing

This reveals a fatal flaw in the authors’ implementation: if the list is infinite, then the recursion will never terminate. For example, (-1) !! [0..] will never terminate, because when the initial index is negative, decrementing the index will never get us to the zero index.

Therefore, I propose the following implementation:

import Prelude hiding ((!!))

(!!) :: [a] -> Integer -> Maybe a
(!!)     [] n = Nothing
(!!) (x:xs) n | n > 0  = (!!) xs (n-1)
              | n == 0 = Just x
              | n < 0  = Nothing

where the first argument is now a list, and the second argument an integer. Note that I used indices of type Integer (“mathematical integers”), instead of type Int (“computer integers”). Let’s see if this implementation works:

*Main> let xs = [7,8,9]
*Main> xs !! 0
Just 7
*Main> xs !! 1
Just 8
*Main> xs !! 2
Just 9
*Main> xs !! (-1)
Nothing
*Main> xs !! 3
Nothing

It appears to be working. No errors. No exceptions. If you, dear reader, happen to be acquainted with Haskell you will almost certainly be shocked (!!!), for this function is trivial! Well, that is true, but I allow myself to be intrigued by trivialities. Moreover, this function is simple enough to allow us to use equational reasoning. For example, let’s compute xs !! 2 using equational reasoning:

[7,8,9] !! 2 = [8,9] !! 1 = [9] !! 0 = Just 9

What if the index is too large? Let’s compute xs !! 4 then:

[7,8,9] !! 4 = [8,9] !! 3 = [9] !! 2 = [] !! 1 = Nothing

Step by step, by successively removing the head of the list, we get where we want to. Unfortunately, this suggests that accessing an arbitrary element of the list will not be \mathcal{O} (1).

__________

References

[1] John O’Donnell, Cordelia Hall, Rex Page, Discrete Mathematics using a Computer (2nd edition), Springer, 2006.

Perfect shuffles in Haskell

April 17, 2012

Consider a deck of 2 n cards. There exists a bijection from the set of cards to the finite set \{1, 2, \dots, 2 n\}, which allows us to represent the deck using a (2 n)-tuple, as follows

(1, 2, \dots, n, n+1, n+2, \dots, 2 n)

where the first element represents the card at the top of the deck, and the (2 n)-th element represents the card at the bottom of the deck. A perfect shuffle consists of cutting the deck exactly in half and then interlacing / interleaving the two halves perfectly. The in shuffle moves the original top card to the second position. The out shuffle leaves the original top card on top [1]. After a perfect in shuffle the deck will be

(n+1, 1, n+2, 2, \dots, 2 n, n),

whereas after a perfect out shuffle the deck will be

(1, n+1, 2, n+2, \dots, n, 2 n).

Is it possible to restore the deck to its original order after a certain number of perfect shuffles? Amazingly (or perhaps not), it is indeed possible!

__________

Example 1

Suppose we are given a deck of 8 cards (i.e., n = 4). A perfect out shuffle of this deck is depicted below

[ source ]

Using pencil and paper, it is easy to see that after three perfect out shuffles the deck will be restored to its original order. Here is a Haskell script:

-- create Card data type
data Card = Ace | Two | Three | Four | 
            Five | Six | Seven | Eight 
            deriving (Show, Eq)

-- create type synonym
type Deck = [Card]

-- interlace two decks of the same size
interlace :: Deck -> Deck -> Deck
interlace [] [] = []
interlace (c1:d1) (c2:d2) = [c1,c2] ++ interlace d1 d2

-- perfect out shuffle
outShuffle :: Deck -> Deck
outShuffle deck = interlace deck1 deck2
                  where deck1 = fst (splitAt 4 deck)
                        deck2 = snd (splitAt 4 deck)

-- perfect in shuffle
inShuffle :: Deck -> Deck
inShuffle deck = interlace deck2 deck1
                 where deck1 = fst (splitAt 4 deck)
                       deck2 = snd (splitAt 4 deck)

where I used function splitAt to cut the deck into two halves of equal size. Since I am not acquainted with any function that performs interleaving, I did create my own function. We load this script into GHCi. Here is a brief GHCi session:

*Main> -- create deck
*Main> let deck = [Ace,Two,Three,Four,Five,Six,Seven,Eight]
*Main> -- perform 3 out shuffles
*Main> take 4 (iterate outShuffle deck)
[[Ace,Two,Three,Four,Five,Six,Seven,Eight],
 [Ace,Five,Two,Six,Three,Seven,Four,Eight],
 [Ace,Three,Five,Seven,Two,Four,Six,Eight],
 [Ace,Two,Three,Four,Five,Six,Seven,Eight]]
*Main>

Indeed, after three perfect out shuffles, we do obtain a deck in the original order! One is tempted to wonder what the order of the deck will be after three perfect in shuffles. Let us see:

*Main> -- create deck
*Main> let deck = [Ace,Two,Three,Four,Five,Six,Seven,Eight]
*Main> -- perform 3 in shuffles
*Main> take 4 (iterate inShuffle deck)
[[Ace,Two,Three,Four,Five,Six,Seven,Eight],
 [Five,Ace,Six,Two,Seven,Three,Eight,Four],
 [Seven,Five,Three,Ace,Eight,Six,Four,Two],
 [Eight,Seven,Six,Five,Four,Three,Two,Ace]]

After three perfect in shuffles we do not obtain the original deck. Instead, we obtain a reversal of the original deck.

__________

Example 2

Let us now work with the standard deck of 52 cards. People say that eight consecutive perfect out shuffles will restore the deck to its original order. Let us check whether this is true or not.

Here is a new Haskell script:

-- create Value data type
data Value = Ace | Two | Three | Four | Five | Six | Seven | Eight | 
             Nine | Ten | Jack | Queen | King deriving (Show, Eq)

-- create Suit data type
data Suit  = Club | Diamond | Heart | Spade deriving (Show, Eq)

-- create type synonyms
type Card = (Value,Suit)
type Deck = [Card]

-- create list of card values
values :: [Value]
values = [Ace,Two,Three,Four,Five,Six,Seven,
          Eight,Nine,Ten,Jack,Queen,King]

-- create list of card suits
suits :: [Suit]
suits = [Club,Diamond,Heart,Spade]

-- create deck
deck :: Deck
deck = [(v,s) | v <- values, s <- suits] 

-- interlace two decks of the same size
interlace :: Deck -> Deck -> Deck
interlace [] [] = []
interlace (c1:d1) (c2:d2) = [c1,c2] ++ interlace d1 d2

-- perfect out shuffle
outShuffle :: Deck -> Deck
outShuffle deck = interlace deck1 deck2
                  where n = div (length deck) 2
                        deck1 = fst (splitAt n deck)
                        deck2 = snd (splitAt n deck)

-- perfect in shuffle
inShuffle :: Deck -> Deck
inShuffle deck = interlace deck2 deck1
                 where n = div (length deck) 2
                       deck1 = fst (splitAt n deck)
                       deck2 = snd (splitAt n deck)

We load this script into GHCi. Here is a GHCi session:

*Main> -- print deck
*Main> deck
[(Ace,Club),(Ace,Diamond),(Ace,Heart),(Ace,Spade),
 (Two,Club),(Two,Diamond),(Two,Heart),(Two,Spade),
 (Three,Club),(Three,Diamond),(Three,Heart),(Three,Spade),
 (Four,Club),(Four,Diamond),(Four,Heart),(Four,Spade),
 (Five,Club),(Five,Diamond),(Five,Heart),(Five,Spade),
 (Six,Club),(Six,Diamond),(Six,Heart),(Six,Spade),
 (Seven,Club),(Seven,Diamond),(Seven,Heart),(Seven,Spade),
 (Eight,Club),(Eight,Diamond),(Eight,Heart),(Eight,Spade),
 (Nine,Club),(Nine,Diamond),(Nine,Heart),(Nine,Spade),
 (Ten,Club),(Ten,Diamond),(Ten,Heart),(Ten,Spade),
 (Jack,Club),(Jack,Diamond),(Jack,Heart),(Jack,Spade),
 (Queen,Club),(Queen,Diamond),(Queen,Heart),(Queen,Spade),
 (King,Club),(King,Diamond),(King,Heart),(King,Spade)]
*Main> -- perform 8 out shuffles
*Main> let deck8 = (iterate outShuffle deck) !! 8
*Main> -- is the new deck the same as the original one?
*Main> deck8 == deck
True

Indeed, it is true! Eight perfect out shuffles suffice!

__________

References

[1] Persi Diaconis, R. L. Graham, William M. Kantor, The Mathematics of Perfect Shuffles, Advances in Applied Mathematics, Volume 4, Issue 2, June 1983.

Feedback systems in Haskell

January 30, 2012

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.

__________

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 \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!

Cascading systems in Haskell

January 25, 2012

Last weekend we learned how to build systems (LTI or otherwise) using Haskell. We now want to construct interconnections of systems. In this post we will study the series interconnection of systems, usually known as cascade interconnection. Parallel and feedback interconnections will be discussed in future posts.

__________

Cascading two LTI systems

Let us consider the series interconnection of two causal discrete-time LTI systems, \mathcal{H}_1 and \mathcal{H}_2, as depicted below

where y = \mathcal{H}_1 (x) and w =\mathcal{H}_2 (y) are the outputs of each LTI system in the cascade. Since the output of system \mathcal{H}_1 is the input of system \mathcal{H}_2, we have w = (\mathcal{H}_2 \circ \mathcal{H}_1) (x).

What LTI systems should we consider? Let us choose the simplest ones: the accumulator and the differentiator. Thus, let \mathcal{H}_1 be an accumulator (also known as “discrete-time integrator”), whose input-output relationship is as follows

y (n) = y (n-1) + x (n),

and let \mathcal{H}_2 be a first difference operator (also known as “discrete-time differentiator”), whose input-output relationship is

w (n) = y (n) - y (n-1).

Note that the output of the cascade of these two LTI systems is thus

w (n) = y (n) - y (n-1) = (y (n-1) + x (n)) - y (n-1) = x (n)

and, hence, the cascade is input-output equivalent to the identity operator. Since (\mathcal{H}_2 \circ \mathcal{H}_1) (x) = x for all signals x, we say that \mathcal{H}_2 is the left-inverse of system \mathcal{H}_1. Since both systems are LTI, the operators commute, i.e., \mathcal{H}_2 \circ \mathcal{H}_1 = \mathcal{H}_1 \circ \mathcal{H}_2 and, therefore, \mathcal{H}_2 is also the right-inverse of system \mathcal{H}_1. Since \mathcal{H}_2 is the left- and right-inverse of \mathcal{H}_1, we say that \mathcal{H}_2 is the inverse of system \mathcal{H}_1 (and vice-versa). Do keep in mind, however, that not all systems are invertible. For details, take a look at Oppenheim & Willsky [1].

__________

Implementation in Haskell

We can easily find a state-space realization for the accumulator, but that will not be necessary. As in previous posts, let us view discrete-signals as lists. Thus, the accumulator takes a list [x_0, x_1, x_2, \dots], and returns the following list

[y_0, y_1, y_2, \dots] = [x_0, x_0 + x_1, x_0 + x_1 + x_2, \dots]

where we assume that the initial condition of the accumulator is zero (i.e., y_{-1} = 0). Instead of using scanl yet once again (which would require us to drop the head of the list), let us now use scanl1 to implement the accumulator:

acc :: Num a => System a a
acc = scanl1 (+)

Please note that if the initial condition of the accumulator is not zero, we should use the following code instead:

acc' :: Num a => System a a
acc' us = tail $ scanl (+) acc_ini us

where acc_ini is the initial condition of the accumulator (analogous to the constant of integration in integral calculus).

The differentiator is not a proper system [2] and, therefore, it has no state-space realization. The differentiator takes a list [y_0, y_1, y_2, \dots], and returns the following list

[w_0, w_1, w_2, \dots] = [y_0, y_1 - y_0, y_2 - y_1, \dots]

where we again assume that y_{-1} = 0. Note the following

[w_0, w_1, w_2, \dots] = [y_0, y_1, y_2, \dots] - [0, y_0, y_1, \dots]

i.e., list w is obtained by elementwise subtraction of a right-shifted version of list y from list y itself. Subtracting two lists elementwise can be implemented using zipWith:

diff :: Num a => System a a
diff ys = zipWith (-) ys (0 : ys)

If the initial condition of the differentiator is not zero, we should instead use the following code:

diff' :: Num a => System a a
diff' ys = zipWith (-) ys (diff_ini : ys)

where diff_ini is the initial condition of the differentiator. Piecing it all together, we finally obtain the following Haskell script:

type Signal a = [a]
type System a b = Signal a -> Signal b

-- accumulator
acc :: Num a => System a a
acc = scanl1 (+)

-- differentiator
diff :: Num a => System a a
diff ys = zipWith (-) ys (0 : ys)

-- cascade of the acc. and diff.
sys :: Num a => System a a
sys = diff . acc

Take a look at the last line. It says that cascading systems is the same as composing systems! Hence, the Haskell implementation is conceptually very close to the mathematical formulation using operators. Functional analysis meets functional programming

We run the script above on GHCi and then play with it:

*Main> -- build unit impulse
*Main> let delta = 1.0 : repeat 0.0 :: Signal Float
*Main> -- output of the accumulator
*Main> let ys = acc delta
*Main> take 20 ys
[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> -- output of the differentiator
*Main> let ws = diff ys
*Main> take 20 ws
[1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
*Main> -- impulse response of the cascade
*Main> let hs = sys delta
*Main> take 20 hs
[1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

The impulse response of the accumulator is the unit step. The impulse response of the cascade is the unit impulse, as we expected. The differentiator is the inverse of the accumulator, and vice-versa.

__________

References

[1] Alan V. Oppenheim, Alan S. Willsky, S. Hamid Nawab, Signals & Systems, 2nd edition, Prentice-Hall, 1997.

[2] Panos Antsaklis, Anthony Michel, A Linear Systems Primer, Birkhäuser Boston, 2007.


Follow

Get every new post delivered to your Inbox.

Join 77 other followers