Animation with MATLAB: Random walk in 2D space

This tutorial will show you how you can make a small animation to show the random walk in 2D space. Like other programs, animation is really easy. Let us see how we can do it.

Random Walk

Random walk is a walk or series of steps where each steps taken by the object is independent of previous steps. It is also known as drunken walk where the person does not know where to go and each steps taken by him is random.

We will consider the discreet case in 2D space and the steps allowed are only four as shown in the figure 1.

Fig 1 : 2D space: Possible directions for random walk

At each time step, the object chooses a number randomly out of {1,2,3,4} and depending on that number, it will choose one of the four direction and take a step. Similar process is followed at each time step and this continues for infinite time. Random walk is the basis of many natural process like Brownian motion.

Simulation of Random Walk

Since it is an iterative process, we will need to use a loop for simulation of this. We will first simulate a single step then we will put a loop around it.

So we need to put four things.
1. assign the current position r = [x,y]

r=[0 0];

2. Get a random number B to select the direction
We will generate a number B with two binary values [b1, b2].
This two values combined will select the direction as following
\[ \begin{array}{|r|c|l|L|}
B & & Direction &  rStep \\
----- & ----- & -------------- & ---------------- \\
[1 ,1] &-> &+X direction & [1 ,0] \\ 
[1 ,0 ] &-> &+Y direction & [0 ,1]\\
[0 ,1 ] &-> &-X direction &[-1 ,0]\\
[0, 0 ] &-> &-Y direction &[0 ,-1]\\
\end{array}
\]

B1=rand(1,2);
B=B1>0.5;

B will be 0 or 1 depending on whether B1 is less than 0.5 or not.

3. calculate new position of object with the step in the selected direction.
    rnew = r + rStep;

    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end

Remember if you put  B==1 as condition inside if, it will be considered true only when all the elements in B satisfies the condition.

So rnew=r+[1 0] only when B(1) and B(2) both are 1. 

Otherwise if one of them is 0, it will check which one is 1 and decide on. In the end, if no conditions are met, it means that B is [0 0] and negative Y direction is chosen.

4. Update the position with this new value

%UPDATE the new position
r=rnew;

Now let us put together

    r=[0 0];
    B1=rand(1,2);
    B=B1>0.5;
    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end
    %UPDATE the new position
     r=rnew;

Now put a Loop around it to make it an iterative process.

r=[0 0];
for t= 0:0.1:100
    B1=rand(1,2);
    B=B1>0.5;
    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end
   
    %UPDATE the new position
     r=rnew;
end

Here the step size is 1 and time step is 0.1 so the speed of the walk is 10 unitdistance/unittime

Animation of the walk

Only thing left is now to show the walk at each time. We will plot the step by a line which will connect the previous r to the current rnew.

plot([r(1) rnew(1)],[r(2) rnew(2)]);

Since you are adding each step on the previous trace, you need to use hold on atleast once so that previous steps are not deleted.

hold on;
plot([r(1) rnew(1)],[r(2) rnew(2)]);

Now the complete code looks like t a Loop around it to make it an iterative process.

r=[0 0];
for t= 0:0.1:100
    B1=rand(1,2);
    B=B1>0.5;
    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end
   
    hold on;
    plot([r(1) rnew(1)],[r(2) rnew(2)]);
    %UPDATE the new position
     r=rnew;
end

when you run the above code, it does not give any output. The missing thing is the command 'drawnow'

Draw Now

When MATLAB sees a series of equations(commands) and a plot command in between, it calculates the graphics but it ignores the plot command at the moment to save time and when all the computation is done, then only it shows the plots. But since you want the steps to be shown at each time, we can force MATLAB to flush the plot at the moment by using drawnow command.

So the final code looks like as

r=[0 0];
for t= 0:0.1:100
    B1=rand(1,2);
    B=B1>0.5;
    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end
   
    hold on;
    plot([r(1) rnew(1)],[r(2) rnew(2)]);
    drawnow ;
    %UPDATE the new position
     r=rnew;
end

Run this code and you will get the following plot figure 2, You can increase time limit by changing the for loop time vector.

Figure 2 : Random Walk

Adding Drift

In the plot you can see that, as time passes object starts coming back to origin. It never stays at origin but it has the tendency to move around the origin and Expected translation at t= infinity is 0. This is because selection of each direction is eqaully probable since you compared with 0.5. So

\[

Pr[Selected Direction = +X ] = P[B(1)=1 \& B(2)=1] \] \[=P(b_1 =1) \times P(b_2 =1) = P(b>0.5) \times P(b >0.5)=\frac{1}{4}

\]

Similarly is the other cases.

Now let us change the probabilities and see the effect.

r=[0 0];
for t= 0:0.1:100
    B1=rand(1,2);
    B=B1>[0.4 0.3];
    if B==1
        rnew=r+[1 0];
    elseif B(1)==1
            rnew=r+[0 1];
    elseif B(2)==1
        rnew=r+[-1 0];
    else
        rnew=r+[0 -1];
    end
   
    hold on;
    plot([r(1) rnew(1)],[r(2) rnew(2)]);
    drawnow ;
    %UPDATE the new position
     r=rnew ;
end

Now we are comparing b1 with 0.4 and b2 = 0.3. So the probability becomes

\[

Pr[Selected Direction = +X ] = P[B(1)=1 \& B(2)=1]\] \[

=P(b_1 =1) * P(b_2 =1) = P(b>0.4)*P(b >0.3)=0.6*0.7=0.42

\]

Similarly with other cases.

The result is shown in figure 3. You can clearly observe a tendency to move in +XY directions. This shows the motion of the object when a drifting force exists in the medium.

Fig 3: Random walk with drift

Simulation of a closed loop system with a Controller

In the previous post http://matlabbyexamples.blogspot.com/2011/07/simulation-of-system-described-by-its.html, we have seen the response of the system was unbounded even for bounded input (a step signal). Here we will try to put a controller to stablize the output.

A feedback with controller

We will use a standard technique to stabilize the output known as feedback where output of the plant is fedback to the plant as an input. The Figure 1 illustrates an system and a derived system formed by connecting output to input through a controller. Controller is a small system which takes the difference between original input and the output and generates a controlled input which is fed to main system. This derived system is also known as closed loop system while the original is known as open loop system.

Fig 1 : Open and Closed Loop System

Building an Environment 

First we will make a controller which is just proportional controller with gain 10.

function v = controller1(e)
k=10;
v=k*e;

 Now build an modified environment which contains the closed loop model and its stimulus input. Let us take step signal as input for this example.

function dybydt = env2(t,y)
%Reference Input Signal
u=1*(t>0);
%Input to controller: Feedback from sys1
e=u-y;
%Controller
v=controller1(e);
%Plant

dybydt = sys1(t,y,v) ;

Simulating the system

The above build environment can be simulated using ode45. Create a blank script file and write the following in that.

%define timespan
tspan=[0 10];
%define initial value
y0=8;
[t y]= ode45(@env2,tspan,y0);

 The output is defined as y and can be plotted by

plot(t,y);

as in Fig 2.

Fig 2: Response of a Closed Loop System

 Now you can see, output is stabilized and bounded for step input. You can explore the effect of controller by changing the value of controller gain k and type of controller.

Simulation of a system described by its Dynamics

The present post and following few posts will talk about the simulation of any system described by its dynamics in MATLAB. We will see how we can simulate an open system and closed system with some small examples.

System Dynamics

Any system is described by its dynamics, Dynamics, in simple words, means how the outputs are dependent on current input and some internal variables called state along with how derivatives of  these states are behaving. Take this example of an electronic circuit with a Resistance R, an inductance L and a capacitance C connected as in Figure 1. Output is the voltage across L and input is the voltage source. This system has two internal parameters Vc and i which affect the output as

 \[ V_o = V-(V_c + iR) \]

Now this Vc and i are dynamic and their derivatives are given as

\[\frac{dV_c}{dt}=\frac{i}{C}\]

\[\frac{di}{dt}=\frac{1}{L}(V-(V_c + iR))\]

As you can see derivatives of states are dependent on state and inputs and so is output. This type of representation is also called state space model.

Fig 1: LCR System with two states

For simulation purposes, we will take two models- one is simple integrator system and second is cruise system.

System1-Integrator model

This system has only one state y which is its output too. State dynamics are given as

\[ \frac{dy}{dt}=u(t)\]

where u(t) is the input to this system.

To create this model, we will make a function sys1 with input u, y and t and which returns derivative term of y. This function describes the dynamics of system, or in other words it represents the system itself.

function dybydt = sys1(t,y,u)

dybydt = u ;

Building an Environment

 Now build an environment which contains the model and its stimulus input. Let us take step signal as input for this example.

function dybydt = env1(t,y)
u=1*(t>0);

dybydt = sys1(t,y,u) ;

Simulating the system

The above build environment can be simulated using ode45. Create a blank script file and write the following in that.

%define timespan
tspan=[0 10];
%define initial value
y0=8;
[t y]= ode45(@env1,tspan,y0);

 The output is defined as y and can be plotted by

plot(t,y);

as in Fig 2.

 

Fig 2: Output of system 1

Understanding the Effect of Sampling on Fourier Domain

This solution is answer to an excercise 7.23 of the Book Openheim, Signal and Systems.

Fig1 : Block Diagram of System

As we can see from the diagram
\[ x_p(t) = x(t) \times p(t) \]
in Fourier domain, it can be written as
\[ X_p(j\omega) = \frac{1}{2\pi}X(j\omega) * P(j\omega) \]
Now calculation the fourier transform of p(t) is little tricky. It is a periodic signal with period 2T. Let us write it as the same form as section 7.1.1
\[ p(t) = \sum_{n=-\infty}^{n=\infty}{(-1)^n \delta(t-n\Delta)} \]
Its fourier transform is given as
\[P(j\omega)= \sum_{n=-\infty}^{n=\infty}{(-1)^n e^{(-n\Delta\omega)}}\]
Now we will simplify it using the fact that (Eq 1)
\[ \sum_{n=-\infty}^{n=\infty}{e^{(-n\Delta\omega)}}= \frac{2\pi}{\Delta} \sum_{n=-\infty}^{n=\infty}{\delta\left(\omega-n\frac{2\pi}{\Delta}\right)}\]
Now ,
 \[P(j\omega)= \sum_{n=-\infty}^{n=\infty}{(-1)^n e^{(-n\Delta\omega)}}\] \[P(j\omega)= \sum_{m=-\infty}^{m=\infty}{(-1)^{2m} e^{(-2m\Delta\omega)}+(-1)^{2m+1} e^{(-2(m+1)\Delta\omega)}}\]  \[P(j\omega)= \sum_{m=-\infty}^{m=\infty}{e^{(-m(2\Delta)\omega)}- e^{(-1\Delta\omega)}e^{(-m(2\Delta)\omega)}}\]\[P(j\omega)= \sum_{m=-\infty}^{m=\infty}{e^{(-m(2\Delta)\omega)}}-e^{(-1\Delta\omega)} \sum_{m=-\infty}^{m=\infty}{  e^{(-m(2\Delta)\omega)}}\]
using the fact above stated (Eq 1)
\[P(j\omega)=\frac{2\pi}{2\Delta} \sum_{m=-\infty}^{m=\infty}{\delta\left(\omega-m\frac{2\pi}{2\Delta}\right)}-e^{(-1\Delta\omega)} \sum_{m=-\infty}^{m=\infty}{  \delta\left(\omega-m\frac{2\pi}{2\Delta}\right)}\]\[P(j\omega)= \frac{2\pi}{2\Delta} \left( 1-e^{(-1\Delta\omega)}\right) \sum_{m=-\infty}^{m=\infty}{\delta\left(\omega-m\frac{2\pi}{2\Delta}\right)}\].
 So it is basically scaled version of impulse train.

Now remaining is to convolve this P(jw) with X(jw) (Refer section 7.1.1). Which results in repitition of X(jw) at 1/2delta interval.

Fig 2: Effect of sampling on Fourier Transform

Fitting Autoregressive Model into the Experimental/Plant Data

In this tutorial, we will learn how we can fit an autoregress model to an experimental data. For our case we generate data from a sample plant with some transfer function and fit a first order ARMAX model to it and compare the results.
Experimental Plant
We will use a plant with transfer function
\[G(z)=\frac{0.333}{z+0.6666}\]
Write following code in MATLAB to generate the model.

g=tf(1/3,[1 .66],0.1)
disp(g);

MATLAB will make a discreet model with name g and display it as

Transfer function:
 0.3333
--------
z + 0.66

Sampling time: 0.1

Simulation of the Plant
We will feed step input to this plant and add some noise to the output to make it more real.
Write following code into MATLAB

%simulate the system
[y t]=step(g,10); 

%add some random noise
n=0.001*randn(size(t));
y1=y+n;

%plot the data with noisy datastairs(t,y);
hold on
stairs(t,y1,'r');

x=ones(size(t)); %input step

ARMAX model
Let us fit a first order model to this data-set [x,y] with the following form

y[n+1]=a y[n] + bx[n]

This equation can be modified as
\[ \left[ \begin{array}{c}y[2]\\y[3]\\ ... \\ y[n+1] \end{array}\right]=  \left[ \begin{array}{cc}y[1] & x[1] \\y[2] & x[2]\\ ... & ... \\ y[n] & x[n] \end{array}\right] \times \left[ \begin{array}{c}a \\b \end{array}\right]\]
or

Y[n+1]=[Y[n]   X[n]] x [a b]'

Y[n+1]=R[n] x Theta

This can be solved for Theta=[a b]'  easily using least square fit method by

Theta= Inv(Y[n+1])*R[n]

So write in MATLAB the following script

Rmat=[y1(1:end-1) x(1:end-1)];
Yout=[y1(2:end)];

P=Rmat\Yout;
disp(P);

You will see the Parameters P as

 -0.6606
  0.3336  

so the fitted model is

 
y[n+1]=-0.6606 y[n] + 0.3336 x[n]

Error:

Let us calculate the fitting error in the model. The error comes due to inserted noise in the output. In real system, this can be anything: process noise, observation noise. So error is nothing but the difference between actual output and output from the fit model using the calculated value of a and b.

%simulated output
Yout1=Rmat*P;

%Mean Square Error
e=sqrt(mean((Yout-Yout1).^2));
disp(e);

This error e is the root mean square error (RMSE) which comes out to be

0.0013572                   

 which is close to the variation of the random noise we add added.

Converting a continuous statespace model to ARMA model

In many situation we need to convert a continuous statespace model in to a linear difference equation model or ARMA/ARMAX discreet time model. This post describes this conversion with some simple examples. 

Example 1: Conversion of First Order Differential Equation to Linear Difference Equation


we will start with a simple example here.  
Problem:
Let us take a simple first order differential equation: 

   dx
  -----  =  - x
   dt

This state space model is with A=-1 and B=0,

we want to discreetize it with T=1s and get a linear differential equation model;

Solution:

First calculate A1 and B1 as following:

1. Using the properties of z=e^-(st), we get that

A1=exp(AT);

A1=1/e;

 or in other words, we can say that

  A1=Linv((sI-A)^-1)  with t=T; 
  A1=Linv(1/(s+1))     =exp(-t) at (t=T)   =e^-1

then similarly

  B1= A^-1 (A1-I)B;= -1(e-1)*0 = 0

 

so the new model is

x[n+1]=A1*x[n]+B1*u[n]

x[n+1]=1/e x[n]   

So the final model in linear differential equation is

  
           x[n+1]-1/e x[n]=0

is the ARM model you wanted

Note :
       exp used is matrix exponential
       exp(A)= I+A+A^2/2!+A^3/3!+... = Linv{(sI-A)^-1}

Example 2: Considering the Inputs

Problem:
Now consider this

dx               
---- =-x+u(t);
dt               

Solution:

So we have

 A1=exp(AT);

 A1=1/e;

 or in other words, we can say that

  Linv((sI-A)^-1)  with t=T; 
 A1=Linv(1/(s+1))     =exp(-t) at (t=T)   = e^-1

then

  B1= inv(A)(A1-I)B;=-1(1/e-1)*1=1-1/e;

so the new model is

x[n+1]=A1 x[n]+B1u[n]

       x[n+1]=1/e x[n]+(1-1/e)u[n];

        x[n+1]-1/e x[n]=(1-1/e)u[n];

x[n+1]-1/e x[n]=(1-1/e)u[n];

Another Method Using Transfer Functions

Let us learn one more method for this.

first calculate the transfer function of this state space model
which is

 G(s)= 1/(s+1);

now discreetize it by using

z=exp(sT)=exp(s) as T=1;

This pretty tricky so we will take the help of MATLAB here to calculate it.

On MATLAB, write the following :

 %make the transfer function using tf
 g=tf(1,[1 1]);

 %convert it to discreet
 g1=c2d(g,1);

you will get

g1=(1-1/e)/(z-1/e);

Now it is very easy to convert it to time domain

X(z)/U(z)=(1-1/e)/(z-1/e);

z X(z)-X(z)/e=(1-1/e)U(z);

The final model is

x[n+1]-1/e x[n]= (1-1/e)u[n]

which is same as earlier one.

Making a MATLAB Media Player

Recently i got a request to write a blog on making a media player on MATLAB. This example shows you a simple way to make a mp3 player in steps. In this step, i will write how to setup the initial background and design base GUI using GUIDE.


Using a mp3 reader
MATLAB does not have an inbuilt mp3 reader function. mp3 is a  encoded which you need to convert to wavfile before reading it. But Dont worry, matlab central provides you all.

go to this download link and download the zip file.
http://www.mathworks.com/matlabcentral/fileexchange/6152-mp3write-and-mp3read

extract the folder and save it to your current directory with name mp3readwrite.

Now you need to add the path of these files in your matlab file.

path(path,'mp3readwrite');

To verify that you have done right, put a mp3 file with name playfile.mp3 in your current directory and run this command

[y f]=mp3read('playfile.mp3');

if you get undefined function or variable mp3readwrite, then check the path of the file inside folder and change the path accordingly.

Generating gui fig using GUIDE
first step is to generate a fig file which contains the basic component of a media player. (As in fig 1)

Fig 1 : sample media player GUI

1. A playlist (ListBox)
2. An Add File button (PushButton)
3. A Volume Silder (Slider)
4. A play Progress Slider(Slider)
5. Play/Pause Button (Push Button)
6. Stop Button (Push Button)
7. Visualization (Axis)

See fig 2 to see the sample tag name used in this example
save this file as mmp.fig and this will generate a mmp.m file

Fig 2 : GUIDE TAGS

 

 In the previous section, We have seen the basic GUI fig. In this section, we would be learning to add basic features in this player.

Making the playlist work
Making the playlist is easy task. Since we need the full path for playing purpose and we need to show only title information in playlist. So we will maintain a spearate list playfiles for full path.

Add this to your opening function 

function mmp_OpeningFcn(hObject, eventdata, handles, varargin)

handles.playfiles=[];

% Update handles structure

guidata(hObject, handles);

We have made the add button. Modify the callback for pladd 

function pladd_Callback(hObject, eventdata, handles)

% hObject    handle to pladd (see GCBO)

% eventdata  reserved - to be defined in a future version of MATLAB

% handles    structure with handles and user data (see GUIDATA)

[f g]=uigetfile;

p=get(handles.playlist,'String');

p{length(p)+1}=[f];

set(handles.playlist,'String',p);

handles.playfiles{length(handles.playfiles)+1}=[g '\' f];

guidata(hObject,handles);

As you can see, uigetfile is used to open a filebrowser and then i store the filename and its path into a list playfiles.
I also update the playlist string property with the file name which is shown on playlist on player.

How to play Audio
There are some option like wavplay or sound in matlab which can play your audio. But you cannot pause or resume it once you have started it. We will use audioplayer class of MATLAB for our purpose.

This small code tells you how to play a wav file by audioplayer object.
first read the .wav file using wavread

[y f]=wavread('aa.wav');

now make an instance of audioplayer with this y and f.

pl=audioplayer(y,f);

pl is the handle to this audioplayer.
Now you can use the play/pause/resume/stop methods to do the respective operations. As an Example:

play(pl);
pause(pl);

resume(pl);

You can get the total sample and current sample number using get method over pl object as 

cs=get(pl,'CurrentSample');

ps=get(pk,'TotalSamples');

To get the progress of the sound file, you can use

ctime=cs/ps;

you can use isplaying option to check whether the pl is running or not.

f=isplaying(pl);

You can combine this with mp3reader to read a mp3 file

cfile='aa.mp3';

extn=cfile[end-2:end];

switch (extn)

  case 'mp3'

        [y f]=mp3read(cfile);

  case  'wav'

       [y f]=wavread(cfile);

end

pl=audioplayer(y,f);

play(pl);

There is one more cool thing to do here. You can make a function which gives you the current progress of music and it will be called automatically by matlab after a certain time if you give its name to audioplayer object. let us see how we can do it.
first make an new m file and save it as UpdateWin.m and write the following in this

function UpdateWin(pl)
%global pl;
c=get(pl,'CurrentSample');
t=get(pl,'TotalSamples');
disp(c/t);

Now make a function f where you can pass this pl.

f=@() UpdateWin(pl);

Now set this f as timer function of pl.

set(pl,'TimerFcn','f()');

Timercall time is automatically set as TimerPeriod: 0.0500 in seconds


Whole piece of code will look like as

cfile='aa.mp3';
extn=cfile[end-2:end];
switch (extn)

  case 'mp3'
        [y f]=mp3read(cfile);

  case  'wav'
       [y f]=wavread(cfile);
end
pl=audioplayer(y,f);

f=@() UpdateWin(pl);
set(pl,'TimerFcn','f()');
play(pl);

Adding Playing Callback
In playbuttonCallback function add the following code for starting the playing of music. Steps we are following are
1. Selecting the file name from playlist

v=get(handles.playlist,'Value');
cfile=handles.playfiles{v};

2. create the player and assign it to the player in handles structure.

extn=cfile[end-2:end];
switch (extn)
  case 'mp3'
        [y f]=mp3read(cfile);
  case  'wav'
       [y f]=wavread(cfile);
end
pl=audioplayer(y,f);
f=@() UpdateWin(pl);
set(pl,'TimerFcn','f()');
handles.pl=pl;

3. Now play this if the player is not running else pause it

if isrunning(pl)
    pause(pl);
else
  resume(pl);
end

Adding other callbacks
making other call backs are pretty easy.
just write following in stop callback

stop(pl);

Making the playersilder move
 Noe it is time to move your player slider to show the progress.
Update your UpdateWin file as following

function UpdateWin(pl)
c=get(pl,'CurrentSample');t=get(pl,'TotalSamples');
%disp(c/t);
h=guidata(pl);
set(h.playslider,'value',c/t);

Enjoy Your player is working now