© Janne Koljonen

University of Vaasa/School of Technology and Innovation Updated: February 3, 2020

ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Introduction to DSP 2

Summary
This document is an introduction to linear filters and dynamical systems. Filters have later a dedicated
lecture. Dynamical systems are an important tool that is used in control theory and systems modelling.
Instead of always going through rigor mathematical definitions and proofs, many things are introduced in
the practical level and with animations for a quick glance and grasping the basic idea.

Table of contents
1. Introduction..........................................................................................................................................2
2. Convolution and filtering.....................................................................................................................3
2.1. Filtering.......................................................................................................................................5

2.2. FIR filter.....................................................................................................................................6

2.3. Image filters................................................................................................................................7

2.4. IIR filters.....................................................................................................................................9

3. Dynamical systems.............................................................................................................................11
3.1. State space model......................................................................................................................15

3.1.1. General linear state space model.....................................................................................17

3.1.2. State space model simulation in Matlab..........................................................................18

3.1.3. Analysis of linear state space model...............................................................................21

3.2. IIR filter as an SSM..................................................................................................................22

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

1. Introduction
Filters are used to modify signals; usually, to remove unwanted components from it. Filters belong to the
body of knowledge as for signal processing. Dynamical systems are an even broader concept. Namely, we
will see that linear filters can be implemented as a linear dynamical model.

The basic filter theory deals with linear filters (see also Figure 1). A linear system H is such that for
signals x and y:

H (ax +by )=aH ( x )+bH ( y ) .

For nonlinear filters, like the median filter, this equation does not hold. Other nonlinear operations are,
e.g., sliding minimum and maximum.

For more examples of linear and nonlinear systems, see: http://www.dspguide.com/ch5/4.htm.

Figure 1. The output of a linear system is the same, no matter if the signals are scaled (by a and b) and
summed before or after the system.

Example 1:

Let H compute the average of two consecutive samples. Let a = 2 and b = 3. Let x = (1, 3, 5) and
y = (2, 2, 1). Test if the linearity equation holds for these (the edges of the signals are now omitted).

The left-hand side: H (ax +by ) …

ax = (21, 23, 25) = (2, 6, 10).

by = (32, 32, 31) = (6, 6, 3).

ax + by = (2+6, 6+6, 10+3) = (8, 12, 13).

H (ax +by )= ( 8+212 , 12+13

2 )
=( 1 0, 12 .5 )
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

And the right-hand side: aH ( x )+bH ( y ) …

H ( x )= ( 1+32 , 3+52 )= ( 2, 4)
H ( y )= ( 2+22 , 2+12 )=( 2 , 1. 5)
aH(x) = (22, 24) = (4, 8)

bH(y) = (32, 31.5) = (6, 4.5)

aH(x) + bH(y) = (10, 12.5).

The left and right-hand sides are equal. Averaging may obey linearity. 

Linearity makes it possible to use more or less simple mathematics to analyze the properties of the filter.
Important mathematical tools are:

 Convolution: output of the filter in time domain

 Z transform: transfer function, black-box model, of the filter; e.g., stability issues

 Fourier transform: steady-state frequency response of the filter; which frequencies are
attenuated/amplified

 Matrix calculus: so-called state space modeling of linear systems.

2. Convolution and filtering

Convolution is an important mathematical tool as for signal processing. In signal processing, convolution
(denoted by operator ) is used, e.g., to compute the resulting output signal y when an input signal x goes
through a linear, time invariant filter h:

y(t) = x(t)  h(t).

More precisely, the filter h is such that computes a sliding weighted average of the input, and that the
weights of the average do not change over time. The output for digital signals is computed as follows:

N−1
y [ n ] =x [ n ]∗h [ n ] = ∑ x [ n−m ] h [ m ]
m=0 ,

where N is the number of weights (coefficient) in filter h. Here 0 stands for the first sample of the signal.

Example 2:

Let a temperature signal be x[n] = (1, 3, 1, 3, 1, 5, 1, 5, 1). Filter x by h[n] = (0.5, 0.5).

There are two coefficients in filter h; thus N – 1 = 1.

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Using the equation above for the first values of y gives:

1
y [ 0 ] = ∑ x [ 0−m ] h [ m ] =x [ 0−0 ] h [ 0 ] +x [ 0−1 ] h [ 1 ] =
m=0
x [ 0 ] h [ 0 ] +x [−1 ] h [ 1 ]=1⋅0 .5+0⋅0 .5=0 . 5.
For n = 1 the output is:

1
y [ 1 ] = ∑ x [ 1−m ] h [ m ] =x [ 1−0 ] h [ 0 ] +x [ 1−1 ] h [ 1 ] =
m=0
x [ 1 ] h [ 0 ] +x [ 0 ] h [ 1 ]=3⋅0 .5+1⋅0 .5=2.
For n = 2 the output is:

1
y [ 2 ] = ∑ x [ 2−m ] h [ m ] =x [ 2 ] h [ 0 ] + x [ 1 ] h [ 1 ] =1⋅0. 5+3⋅0 .5=2.
m=0

For n = 3 the output is:

1
y [ 3 ] = ∑ x [ 3−m ] h [ m ]=x [ 3 ] h [ 0 ] + x [ 2 ] h [ 1 ] =3⋅0 . 5+1⋅0 . 5=2 .
m=0 

Figure 2 shows the whole output of the filter.

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 2. The result of the convolution in example 2. 

In Moodle, there is a Matlab script file that animates this filter. In order to display it start Matlab,
download and open the m file, highlight the code and press F9.

2.1. Filtering
Filters are used to remove noise from signals. Filters are also used to separate different frequencies from
the signal in order change their mutual volume, balance. The latter is probably familiar from music
players where the amount of bass and treble can be adjusted using an equalizer.

When we are dealing with filters that compute convolution, the eye is on frequencies. With sliding
weighted average one can realize such filters as low-pass, high-pass, band-pass, and band-stop. Moreover,
it is possible to computationally obtain a derivative of the signal, among others.

Let us test the filter from example 2 to a signal that has two sines at different frequencies ( f1 = 5 Hz and
f2 = 400 Hz). The effect of the filter is illustrated by computing the Fourier transform spectrum before and
after filtering.

Figure 3 shows the results. It can be concluded that the filter with h = (0.5, 0.5) attenuates the higher of
the two frequency components. So, the filter is a low-pass filter. The filtered signal looks less noisy,
because usually noise occurs at the high frequencies.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 3. Left column: original signal (red) and the filtered signal (green).
Right column: spectra of the original (red) and the filtered signal (green).

2.2. FIR filter

In signal processing, a sliding weighted average filter with a finite number of coefficients h (more
usually, b) is called a Finite Impulse Response (FIR) filter. FIR computes convolution between its
coefficients and the input x.

N−1
y [ n ] =x [ n ]∗h [ n ] = ∑ x [ n−m ] h [ m ]
m=0 .

Figure 4 gives an FIR example with three coefficients. Each coefficient in this case is 1/3, which implies
that the output is the sliding average of every three input samples:

y [ n ] =x [ n ]⋅(1/3)+x [ n−1 ]⋅(1/3)+x [ n−2 ]⋅(1 /3 ) .

In Figure 4, each box with z–1 means a unit delay (in the sense of the properties of Z transform). Unit
delay is equivalent to a shift register where the data are moved one step forward on every clock edge.
Here these clock edges are counted by variable n.

The output of the FIR in Figure 4 is computed given that the current input x[n] = 5, the previous input
x[n–1] = –2, and the input before that x[n–1] = 2:

y [ n ] =5⋅(1/3)−2⋅(1/3 )+2⋅(1/3 )=1 .7 .

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 4. An FIR filter with three coefficients, b = (1/3, 1/3, 1/3), and the computation of a single output.

FIR filters will be discussed later in detail. However, it is beneficial to check an introductory video now:

https://www.youtube.com/watch?
v=aV7cM9nmv0A&feature=iv&src_vid=mCzoxw4j4uA&annotation_id=annotation_103492

2.3. Image filters

FIR filters, and other too, can be generalized to two-dimensional signals. An image is a typical example
of 2D signals.

Example 3:

Let the image X with 5  5 pixels and the image filter mask/kernel H be, respectively:

0 0 0 1 0 H
0 2 4 4 0
1 3 4 3 0 –1 0 1
1 3 3 3 1 –1 0 1
1 0 0 0 0 –1 0 1

Next we can start computing the filtered image Y by sliding the filter over the image. The first position
when the whole filter fist inside the image is as follows (the bold pixels show the location of H):
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

X Y
Position of H in bold

0 0 0 1 0
0 2 4 4 0 7
1 3 4 3 0
1 3 3 3 1
1 0 0 0 0

The result of convolution at this single position is calculated as follows:

Y(1, 1) = X * H =

0 (–1) + 0 (0) + 0 (+1) +

0 (–1) + 2 (0) + 4 (+1) +

1 (–1) + 3 (0) + 4 (+1) = 7

In image processing, the result of convolution is usually inserted to the center position of the filter mask
(see above). Then the mask is shifted to the next position to obtain a new output pixel:

X Y
Position of H in bold

0 0 0 1 0
0 2 4 4 0 7 3
1 3 4 3 0
1 3 3 3 1
1 0 0 0 0

Here the calculus is:

Y(1, 2) = X * H =

0 (–1) + 0 (0) + 1 (+1) +

2 (–1) + 4 (0) + 4 (+1) +

3 (–1) + 4 (0) + 3 (+1) = 3.

As a last example, the image filter is applied to a real-life image. The results in Figure 5 show that this
image filter picks up vertical edges in the image.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 5. Upper: original image. Lower: filtered image. 

2.4. IIR filters

IIR filter is an extension to FIR filter. Whereas FIR computes weighted sliding average of the inputs, IIR
computes weighted sliding average of both input and outputs. That means that the output is fed back to
the system. This feedback can result in unstable behavior if the gain of the feedback path is too high.
(More precisely: IIR is stable iff all poles of its transfer function are inside the unit circle.)

In general, the output of an IIR filter is:

P Q
y (n )=∑ b i x( n−i)−∑ a j y (n− j )
i=0 j=1 ,

where n is discrete time, x(n–i) and y(n–j) refer to the previous inputs and outputs, respectively. Note that
i starts from 0 (i.e., the current input is also used) while j starts from 1 (i.e., the previous output is the first
term in the feedback sum). a0 is actually assumed to be 1.

Figure 6 depicts an IIR filter with totally five coefficients:

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

2 2
y (n )=∑ b i x( n−i)−∑ a j y (n− j )=b 0 x( n)+b 1 x (n−1)+b 2 x (n−2 )−a 1 y (n−1 )−a 2 y (n−2 )
i=0 j=1
.

Figure 6. IIR filter with three feedforward coefficients bi and two feedback coefficients aj. (By Akilaa
(Own work) [GFDL (http://www.gnu.org/copyleft/fdl.html), CC BY-SA 3.0
(http://creativecommons.org/licenses/by-sa/3.0) or GFDL (http://www.gnu.org/copyleft/fdl.html)], via
Wikimedia Commons) IIR filter is a recursive filter (https://en.wikipedia.org/wiki/Recursive_filter).

Example 4:

Let us say that the coefficients of an IIR filter are b = (1/3, 1/3, 1/3) and a = (1, –0.5, –0.55). Simulate the
output using the following Matlab code:

% Coefficients
b = [1/3, 1/3, 1/3];
a = [1, -0.5, -0.55];
% Unit impulse as the input
x = [1 zeros(1, 49)];
% Output of IIR
y = filter(b, a, x);
figure
subplot(2, 1, 1)
hold on
stem(x)
xlabel('n')
ylabel('input, x')
subplot(2, 1, 2)
hold on
stem(y)
xlabel('n')
ylabel('output, y')

Figure 7 shows the input (a unit impulse) and the output of this IIR filter. This is an example of an
unstable IIR filter.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 7. An unstable IIR. The output to a unit impulse input diverges towards infinity. 

3. Dynamical systems
Dynamical systems are used to describe time dependent behavior of signals and systems. The main
components of dynamical system models include the state vector of the system and an evolution rule that
tells what is the future state given the current state. The evolution rule can be deterministic or stochastic.
Moreover, it is important to define the initial state.

Dynamical systems can be used to model continuous systems, for instance, mechanical systems, quantum
mechanics, and control systems. Now we mainly restrict ourselves to discrete-time systems. It means that
state is changed only at discrete time instances. Moreover, often the state space is also discrete, and the
respective system is called a quantized-data systems.

Example 5:

Let us do a Matlab example using discrete-time system with two continuous state variables:

x1[k] = the balance of your bank account [euros], where k is discrete time (1, 2, 3, …), and

x2[k] = your ‘shopping mood’.

Let us define some evolution rules:

(1) After each time step k, your money increases due to interest and decreases due to some shopping. The
amount of shopping depends on your ‘shopping mood’.

(2) After each time step k, your ‘shopping mood’ increases if your balance is over 1000 euros and
decreases otherwise. Shopping mood is always between 0 and 100 points, inclusively.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

More specific, the rules could be:

(1)
x 1 [ k +1 ] =(1+ r ) x1 [ k ] −sx 2 [ k ] ,

where r is the interest rate multiplier and s is your shopping rate [euros/mood point].

x2 [ k +1]=¿ { x2 [k ]+ dm, if x1 [ k ]>1000 ¿ ¿¿¿

(2) ,

where dm is your rate of change of the shopping mood, and restrict mood to [0, 100]:

x 2 [ k +1 ] =max(0 ,min( x 2 [ k +1 ] , 100))

For example, let the initial values be x1[1] = 1001 (euros), x2[1] = 50. Furthermore, let the interest rate be
r = 0.05 (5%), your shopping rate s = 2 (euros/mood point), and your rate of mood changes dm = 10.

When k = 2: you get 0.051001 euros = 50.05 euros as interests and you will do shopping by 250 euros
= 100 euros. Consequently, x1[2] = 1001 + 50.05 – 100 = 951.05 (euros). Your bank account had over
1000 euros, and therefore your mood is increased by 10 points to 60 points.

This dynamical system can now be modelled in Matlab as follows:

% Constants
r = 0.05;
s = 2;
dm = 10;
% Initialize vectors
N = 300;
x1 = zeros(N, 1);
x2 = zeros(N, 1);
% Initial conditions
x1(1) = 1001;
x2(1) = 50;

% Simulate deterministic discrete-time dynamical system

for k = 1:N
x1(k+1) = (1+r)*x1(k) - s*x2(k);
if x1(k) > 1000
x2(k+1) = x2(k) + dm;
else
x2(k+1) = x2(k) - dm;
end
% Restrict to [0, 100]
x2(k+1) = max(0, min(x2(k+1), 100));
end

% Plot state behaviors separately

figure
subplot(2, 1, 1)
plot(x1)
hold on
title(strcat('r = ', num2str(r), ', s = ', num2str(s), ', dm = ', num2str(dm)))
xlabel('k')
ylabel('balance [eur]')
subplot(2, 1, 2)
plot(x2)
hold on
xlabel('k')
ylabel('mood')

% Plot the state space behavior as vectors

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

figure
plot(x1, x2)
ylim([0 100])
xlabel('balance [eur]')
ylabel('mood')

Using different initial values lead to different behavior. Figure 8 and Figure 10 show two examples with
different parameters but equal initial conditions. It can be seen that the range of the balance is almost the
same in both cases. However, it would be difficult to predict when there will be an especially low plunge.

Figure 8. Simulation of the dynamical system (bank account balance).

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 9. The path of the state space vector.

It is typical to dynamical systems that there are some orbits like in Figure 9.

Figure 10. Another simulation. Shopping rate increased from 2 to 3.

When the interest rate is increased, finally the earned interests exceed the spent money. Consequently, the
balance starts to increase infinitely. However, the system does not necessarily immediately show its
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

instability (Figure 11), because there is random-like variation, even though the system is totally
deterministic!

Figure 11. Simulation with a higher interest rate (10.766%). 

3.1. State space model

The example above showed a nonlinear dynamical system. It was nonlinear as the rule (2) that updated to
shopping mood was not a linear function.

The canonical theory of dynamical systems assumes linear systems. In linear cases, the system equation
can be presented by matrices and thus analyzed by matrix algebra.

In the case of linear continuous-time systems, the dynamics of the system can be written by differential
functions.

Example 6:

Consider an object that is attached to a wall by a spring (Figure 12).

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Mass of the object:

m Position of the object:
y(t)

Spring constant: External force:

k F = u(t)

Viscous friction coefficient: b y=0

Figure 12. Modelling of a moving object attached to a spring.

According to Newton’s second law mass times acceleration is force, i.e., a force causes acceleration:

m ÿ ( t )=∑ F (t) .

Now three forces are present: the time-variant external force u(t), spring force ky(t), and viscous friction
force that is linearly dependent on the speed of the object b ẏ ( t ).

Be careful with the signs (directions) of the forces: Positive direction of y is set to the right as well as
positive external force. In that coordinate frame, if y is positive, spring force should be negative (i.e., to
the left). Similarly, if movement is to the right, viscous friction should cause a force to the left.

Overall, the dynamic model becomes:

m ÿ ( t )=u ( t )−b ẏ ( t )−ky (t).

This is a second-order differential equation. Usually, dynamical systems are modelled using first order
differential equations. The trick is to use a many state variables as the order of the system is; now two. Let
x1(t) and x2(t) be the position and velocity of the object, respectively. Then the dynamic model can be re-
written as a first-order differential equation:

⇔
u ( t ) b x2 ( t ) k x1 ( t )
m ẋ2 ( t )=u ( t ) −b x 2 ( t )−k x 1 ( t )¿ m ẋ 2 ( t )= − −
m m m
This is actually the evolution rule for state variable x2(t). The respective rule for x1(t) is now simply:

ẋ 1 ( t )=x 2 ( t ) ,

because the differential of the position is velocity.

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

These two evolution rules can be written as a set of equations (notice rearranging and added zero terms):

{
ẋ 1 ( t )=0 x 1 ( t )+ x 2 ( t ) +0 u ( t )
−k b 1
ẋ 2 ( t )= x 1 ( t )− x 2 ( t ) + u ( t )
m m m

Finally, this can be written in the matrix form:

[ ][ ][ ] [ ]
0 1 0
ẋ1 ( t ) x1 ( t )
= −k −b + 1 u(t) ,
ẋ2 ( t ) m m x 2 ( t ) m

where x1 and x2 are the state variables position and velocity, respectively.

The output of this system is the position y(t). The output is also presented by matrices:

y (t )= [ 1 0 ]
[ ]
x 1 (t )
x 2 (t )
+0 u(t )
.

3.1.1. General linear state space model

The above moving object example can be generalized. A linear continuous time-invariant system with n
states, p inputs, and q outputs is defined by two matrix-equations:

{[ ][ ][ ]
ẋ 1 (t) x 1 (t) u1 (t)
ẋ 2 (t) x (t) u (t)
=A 2 +B 2
⋮ ⋮ ⋮
ẋ n (t) x n (t) u p (t)

[ ][ ][ ]
y 1 (t) x1 (t ) u 1(t)
y 2 (t) =C x2 (t ) + D u 2(t) .
⋮ ⋮ ⋮
y q (t ) x n (t) u p (t)

A is called the state (or system) matrix.

B is called input the input matrix.

C is called the output matrix.

D is called the feed-through matrix.

In addition, there can be a noise matrix K if the system is not deterministic.

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Explicit discrete time-invariant state space model is defined as:

x (k +1 )= Ax(k )+Bu(k )
y (k +1 )=Cx (k )+Du(k ).
However, the A, B, C, and D matrices are not equivalent to the continuous-time case!

3.1.2. State space model simulation in Matlab

Example 6 (revisited):

Let us continue the moving object example by simulating it in Matlab:

% Constants
k = 10;
b = 0.5;
m = 1;

% System matrix
A = [0 1; -k/m -b/m];
% Input matrix
B = [0;1/m];
% Output matrix
C = [1 0];
% Feedthrough matrix
D = 0;
% Noise coefficient matrix
K = [0;0];

% Ts = 0 indicates continuous time

SSM = idss(A,B,C,D,K,'Ts',0,'NoiseVariance',0);

% External force u:
samples = 300;
u = zeros(samples, 1);
u(1, 1) = 100;
u(2, 1) = 100;
u(40, 1) = -200;
u(41, 1) = -200;
u(100, 1) = 300;
u(101, 1) = 300;
u(140, 1) = -400;
u(141, 1) = -400;

data = iddata([],u, 'InputName', 'External Force');

% Sampling interval
data.ts = 0.1;

figure
plot(data)
xlim([0 samples*data.ts])

figure
% Initial condition: position = 1, velocity = 0
sim(SSM, data, [1;0]);
% Save: y = output, ysd = output deviation, x = states
[y, ysd, x] = sim(SSM, data, [1;0]);
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 13. The external force in the simuation.

Figure 14. The output of the simulation, the position of the moving object.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 15. The state vector trajectories of the moving object. In free ride, the orbit is a spiral. External
force causes the discontinuities.

Figure 16. Animation of the moving object (works in PowerPoint). 

© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

3.1.3. Analysis of linear state space model

The matrices can be used for the analysis of the state space model and, e.g., to obtain the transfer
function.

Controllability means that it is possible – within finite time – to steer the states from any initial condition
to desired values. In the case of the moving object, it would possible to have such an external force
function that a desired location and velocity were reached.

A linear state space model is controllable if and only if:

rank [ B AB A 2 B … A n−1 B ]=n ,

where n is the number of state variables.

Example 6… Is the moving object then controllable?

test = zeros(size(A));
for column=1:size(A, 2)
test(:, column) = A^(column-1)*B;
end
disp(strcat('rank = ', num2str(rank(test)) ))
if rank(test) == size(A, 1)
disp('Is controllable')
else
disp('Is not controllable')
end

The output of the Matlab code shows that it is controllable. But how to obtain a desired state is not an
easy question in the general case. 

A dual of controllability is observability. Observability implies the internal states can be inferred from the
output trajectory.

A linear state space model is observable if and only if:

[ ]
C
CA
rank C A2 =n
⋮
n−1
CA

Example 6…

test = zeros(size(A));
for row=1:size(A, 2)
test(row, :) = C*A^(row-1);
end
disp(strcat('rank = ', num2str(rank(test)) ))
if rank(test) == size(A, 1)
disp('Is observable')
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

else
disp('Is not observable')
end

The output of the Matlab code shows that the moving object state model is observable. 

Stability can be assessed from the eigenvalues of the state matrix A. Actually, those eigenvalues are the
poles of the transfer function.

Example 6… For the moving object, the stability test is:

poles = eig(A)

poles =

-0.2500 + 2.2220i

-0.2500 - 2.2220i

A continuous time-invariant system is stable if the real parts of the poles are negative. Now they are
(–0.25), and thus the moving object system is stable. 

3.2. IIR filter as an SSM

State space model can be used to describe linear filters. The time-domain presentation of an IIR filter is a
difference equation:

y (k )=b 0 u(k )+b1 u (k−1 )+b2 u (k−2 )−a 1 y (k−1 )−a 2 y (k −2) ,

where u(k) is now the input as x is reserved to the state variables. The order of this equation is N = 2. It
can be replaced by a set of N first-order state update equations:

{x1 (k+1)=−a1 x1(k)−a2 x2(k)+u(k) ¿ ¿¿¿

In addition, the output is obtained by:

y (k )=( b1 −b0 a1 ) x 1 ( k )+ ( b2 −b0 a2 ) x 2 (k )+ x 0

.

The same as a state space model is:

{[ ] [ ][ ] [ ]
x1(k+1) −a1 −a2 x1(k) 1
=
x2(k+1) 1 0 x2(k) 0
+ u (k) ¿ ¿¿¿
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

For the derivation and generalization of this see:

https://ccrma.stanford.edu/~jos/fp/Converting_State_Space_Form_Hand.html.

Example 7:

Let the coefficients of an IIR filter be b = (1, 0.5, 0.25) and a = (1, –0.5, 0.1). What is the state space
model of this filter? Is the filter stable?

b = [1, 0.5, 0.25];

a = [1, -0.5, 0.1];

% Convert transfer function to state space model

[A,B,C,D]=tf2ss(b,a)

% Stability can be assessed either from

% the transfer function
disp(strcat('roots = ', num2str(abs(roots(a))') ))
% zero-pole plot
% or the SSM
disp(strcat('system matrix eigenvalues = ', num2str(abs(eig(A))') ))

Outputs:

A = 0.5000 -0.1000
1.0000 0
B = 1
0
C = 1.0000 0.1500
D = 1

poles =0.31623 0.31623

system matrix eigenvalues =0.31623 0.31623

An IIR filter is stable if all its poles are inside the unit circle. Now the absolute values of the poles are <
1, i.e., they are inside the unite circle. The zero-pole plot shows this, too (Figure 17). We will later return
the zero-plots and other filter analysis.
© Janne Koljonen
University of Vaasa/School of Technology and Innovation Updated: February 3, 2020
ICAT3180 Sound Processing web course
Lecture 3: Introduction to DSP 2

Figure 17. The zero-pole plot of the IIR filter. o = zeros,  = poles. 