Median Filtering using MATLAB

In this blog we will learn how we can do a time domain filtering for an image. We will consider median filtering for this example. We will see two methods - first one is the iterative methhod which is time consumingg. Second methods uses matlab's vectorization and performs well.

Example Image
For this blog, we will take a very short image to understand how filtering is being done. Our sample image is shown in Fig1 (7x7).

Fig1 : Image A (magnified)

A =
8 5 8 0 6 8 7
9 9 1 8 1 6 7
1 9 4 9 7 3 1
9 1 9 6 0 9 4
6 9 7 7 2 0 4
0 9 9 7 0 4 6
2 4 6 3 0 3 7

You can create this by the following command

A=[8 5 8 0 6 8 7;9 9 1 8 1 6 7;1 9 4 9 7 3 1;9 1 9 6 0... 9 4; 6 9 7 7 2 0 4;0 9 9 7 0 4 6;2 4 6 3 0 3 7];

In actual cases you will read a image by using imread

A=imread('C:\users\Myname\codes\image1.jpg');
A=rgb2gray(A);

[l m]=size(A);

Kernel
We would be filtering with 3x3 kernel. It means that for each pixel of the output image, we would take a 3x3 neighborhood around the respective pixel in the original image, and assign the median color value in this neighborhood to the pixel of output image.

Output
Output would be a image of 7x7 too. Let us preallocate it with all zeros.

B= zeros(l,m);

Iterative method:
For each i and j, we first extract a 3x3 matrix around this point as shown in figure.

Fig2: Image Kernel Elaboration

K=A(i-1:i+1,j-1:j+1);

This will select a 3x3 matrix containing elements from i-1,i,i+1 row and j-1,j,j+1 column and assign it to K.

Now calculate median of this K and assign it to B(i,j);

B(i,j)=median(K);

We run a loop on i,j = 2 to 6. We cannot start with 1 because 3x3 kernel should fit inside the image to avoid any matrix index out of bound error. In the same way we cannot have i or j = 7. Note that we can take only those i where i-1, i,i+1 all three should come in the range 1 to 7. So with the j.

for i=2:6

  for j=2:6

    K= A(i-1:i+1,j-1:j+1);

    B(i,j)=median(K);

  end

end

 So for example for i=3,j=4, we take the kernel consisting of 2-3-4 row and 4-5-6 column as shown in figure 2. Elements are 1,4,6,7,7,9,9,9,9 and the median is 7 , so the output image's B(3,4)=7.

 This is done for all the indexes.

Padding:
We will first pad the original image to take care of the corner cases. To pad we will create a new matrix with l+2, m+2 dimension and assign the middle region from i=2 to l+1, j=2 to m+1 from the original image.

Y=zeros(l+2,m+2);
Y(2:l+1,2:m+1)=A;

Now remember rth,sth element of Y is equal to ith,jth element of A where i=r-1,j=s-1 for r,s>2.
Now you can go again i=2:l+1, 2:m+1, but remember to assign it as following

B(i-1,j-1)=median(K);

Vectorized Method:
The iterative method is little slow. MATLAB is optimized for vectorized operation. Let us try the same with vectorized way but it is little tricky. As you know, to calculate median, you need to pass MATLAB a matrix.

First we will create 9 shifted version of A (let us say A1, A2,... A9) so that for any i,j,A1(i,j) to A9(i,j) represent the 9 values in 3x3 neighborhood of A(i,j).
since we are considering i=2 to 6, these 9 matrixes will be 5x5 size;

A1=A(1:5,1:5);
A2=A(1:5,2:6);
A3=A(1:5,3:7);


A4=A(2:6,1:5);
A5=A(2:6,2:6);
A6=A(2:6,3:7);


A7=A(3:7,1:5);
A8=A(3:7,2:6);
A9=A(3:7,3:7);

Now we need to find median of 9 values A1(i,j) to A9(i,j) for all i,j. We first combine them into a 3D matrix and calculate the median in its first dimension.

C(1,:,:)=A1;
C(2,:,:)=A2;
C(3,:,:)=A3;
C(4,:,:)=A4;
C(5,:,:)=A5;
C(6,:,:)=A6;
C(7,:,:)=A7;
C(8,:,:)=A8;
C(9,:,:)=A9;
D=median(C,1);
B=squeeze(D);

And Done! Here is you answer B. This is pretty fast specially for large matrixes. There are only 9 operations (equal to kernel size) are happening, while in iterative one, m*n operations are happening, so if m=n=1000, there are 1000000 operations !!

Multi-variate Model Fitting using Taylor Series Method

In the previous example, we have learnt how to calculate a singlevariate model using a taylor series method when no knowledge about form of relation is known between input and output. This example will extend the same using iterative application of singlevariate model fitting.

We will consider a bivariate model and same can be easily generalized for higher dimensions.

Problem Statement:
We consider a system where there are two inputs/parameters (x and y) which we can change and observe the output (z). Our aim is to find a relation between this output and inputs but again we don't know what relation form(whether it is exponential, logarithmic, linear or trigonometric) exists.

Taylor Series Method
see this article for Taylor Series Method http://matlabbyexamples.blogspot.com/2011/12/singlevariate-model-fitting-for.html.

Iterative Taylor Series Method
We will first make one parameter (y) constant (let us say at y0) and change the other one (x) only. Then we try to fit a taylor series model in to it and we get coefficient b_n|at y=y0 such that  \[ f(x) = \sum_n{b_n x^n} \].
Now we will change the parameter y to y1 and do the same fit again to get a different set of coefficients  b_n|at y=y1. Do it for all value of y and each time you get a 1xN vector b containing b_n's. Now stack them such that  we get a matrix B each row represent b for one value of y. Now observe that ith column (let us call it B_i) represents b_i's for all values of y.

Now it is clear that these columns are only dependent on value of y. So again we can fit a taylor series to each of these columns such that  \[ b_n = \sum_m{c_m y^m} \].  This can be done by solving \[B_n = Y*C\]

So we have the final model as
\[ f(x) = \sum_n{\sum_m{c_{m,n} y^m} x^n} \].

MATLAB Implementation:

function C=taylorseriesfit2D(X,Y,Z);
%assumes X and Y are meshgrid format matrix
for i=1:length(y)
z=Z(i,:)';
B(i,:)=taylorseriesfit1D(X(1,:),z);
end

for i=1:size(B,2)
C(:,i)=taylorseriesfit1D(Y(:,1),B(:,i));
end

where taylorseriesfit1D is given as

function B= taylorseriesfit1D(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7];
B=Xmat\z;

Simulation Data generation:

x=[0:0.1:1]';
y=[0:.1:pi/2]';
[X Y ]=meshgrid(x,y);
Z=exp(X).*exp(Y);
C=taylorseriesfit2D(X,Y,Z)

C will come like for N=7.

1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
0.5000 0.5000 0.2500 0.0833 0.0208 0.0043 0.0006 0.0002
0.1669 0.1669 0.0835 0.0278 0.0069 0.0014 0.0002 0.0001
0.0410 0.0410 0.0205 0.0068 0.0017 0.0004 0.0000 0.0000
0.0092 0.0092 0.0046 0.0015 0.0004 0.0001 0.0000 0.0000
0.0008 0.0008 0.0004 0.0001 0.0000 0.0000 0.0000 0.0000
0.0004 0.0004 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000

Final Guesses 
Again we can guess the model by observing these coefficients and plotting them with respect to m,n (Fig1). First thing we can observer that it is symmetric, therefor exihibit the same relations for both the inputs to ouput.

Fig1: Plot C vrs m,n

Again, we cannot guess, even then these matrixes are sufficient for all numerical purposes.

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

In this example we will talk about how we can fit a model/function in a single variate simulation data. In many practical cases, we get some observation data after changing an input variable and our aim is to find the relation between input and output. If we know the relation form atleast (like y=mx+c or y=d*sin(x)^2), then we can find the values of these coeefficients y,c or d by linear/non linear regression. But what if we dont know the relation form, we dont have any information about whether they follow exponential form or logarthim relation. Here is a trick which we are going to discuss in this example.

Taylor Series fit
Let i/o relation be represented by a function z=f(x). We assume here that f(x) is continuous and differentiable, which is generally the case with practical/real time systems. Now each such function can be represented using a taylor's series
\[ f(x) = \sum_n{b_n x^n} \] we will limit the summation to n=N
Since we know the data x and z=f(x), we can calculate the coefficients bi .
Let \[B=\left[\begin{array}{ccccc} b_0 & b_1 & b_2  & ... & b_N \end{array} \right]^T \] and \[ Xmat =\left[ \begin{array}{ccccc} 1 & x^1 & x^2 & ... & x^N \end{array} \right] \]
so \[ z=f(x)= Xmat B \]

This is also true when x is a column vector mx1 and contains m data points and z contains values of f at those m point.
So B can be calculated as \[B=Xmat^{-1} Z\]

function: taylorseriesfit()
Writing the function in matlab is very easy. we just need to follow the steps as explained above.
Since our function need to matrixs x and z, these will be input to the function and B will be the output of the function. We will take N=4 here. Increasing the value may result in accuracy of your fit.

function B= taylorseriesfit(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x^4];
B=Xmat\z; 

Here x and z are column vectors representing data points.

Generation of Simulation Data 
We will generate a sample Simulation data for our example. In real example, you will get this data from some experiment after changing an input parameter(x) and observing the output(z).

x=[0:0.1:1]';

z=exp(x);

Calling the defined function

B=taylorseriesfit(x,z);
disp(B);

You will get B as

B =
1.0000
0.9988
0.5099
0.1400
0.0696

Now your model is
\[ z= x*B = B(1)+B(2)x+B(3)x^2+B(4)x^3+B(5)x^4 \] \[ z= 1+0.99x+0.50 x^2+0.14x^3+.07x^4 \]

Error
we can calculate the modelling error (MMSE) by comparing the actual output and the simulation output.

z_o = Xmat*B;
e=sqrt(sum((z-z_o).^2));
disp(e);

you will get something like

e= 6.9400e-005

which is pretty good.

Final Observation
Now if we observe the coefficients and plot them with respect with n, we may guess the actual function. Fig1 represents such a graph for N=7 case.

Fig1 : b_i plotted versus i

like for this example, it looks like 1/n!, so this is an exponential relations. Well, if we can guess, it is good. But if we cannot, still we have a relation which is sufficient for all numerical purposes.

Circuit Solver using MATLAB Programming

In this blog we will learn to simulate a circuit using MATLAB Programming which is the basic step to create your own circuit solver. The circuit we are going to solve is shown in Fig 1. Here V is 6V.

Fig 1: Circuit

Create NetList and SourceList
You need to first represent this circuit in terms of matrixes so that you can use them in MATLAB. These two matrixs will describe the network- NetList and SourceList
Let us first give a number to all the nodes in the circuit. Two nodes should be separated by atleast an element.
The considered circuit has 4 nodes and let us choose node 4 as the reference/ground node.

Now create a matrix named as NetList with three columns N1, N2 and R. Each row of this matrix will represent a R resistance element between Node N1 and N2.

Now create a matrix named SourceList with two columns N and V. Each row of this matrix will represent a N node with V voltage. Since only two node voltages are know, this will contain two rows.

Fig 2: NetList and Source List

Create A and B matrix
If Total Nodes are N
A is of N x N
Aij=
        Case 1: if i is not a Source Nodes
                       -1/Rij where Rij is the resistance between i and j (if i!=j)
                       Sum of inverse of all resistance connected to i (if i==j)

        Case 2: if i is a source Node
                        Aii=1 and
                        rest Aik=0

If Total Nodes are N
B is of N x 1
Bi=
       Case 1: if i is not a Source Node
                    0
       Case 2: if i is a source Node
                    Voltage of that source Node



Nodal Equations

Using the above nodal equation can be represented as
\[\mathbf{AV=B}\]\[\mathbf{V}=[V_1 V_2 ... V_N]^T\]\[V=A^{-1}B\]

This V represent the Voltage of all the nodes, thus circuit is solved.


MATLAB Code
Initialization of the matrixs

NetList=[1 2 2; 1 3 4; 2 3 5.2; 3 4 6; 2 4 3];
Vnod=[1 6; 4 2];
l=size(NetList,1);
N=max([NetList(:,1) ;
NetList(:,2)]);
A=zeros(N,N);
 B=zeros(N,1);

A Matrix Creation

for i=1:l
       n1=NetList(i,1);
       n2=NetList(i,2);
       if n1==n2
       else
          A(n1,n2)=A(n1,n2)-1/NetList(i,3);
          A(n2,n1)=A(n2,n1)-1/NetList(i,3);
          A(n1,n1)=A(n1,n1)+1/NetList(i,3);
          A(n2,n2)=A(n2,n2)+1/NetList(i,3);
      end
 end

B Matrix Creation

for i=1:size(Vnod,1)

        A(Vnod(i,1),:)=zeros(1,N);

        A(Vnod(i,1),Vnod(i,1))=1;

        B(Vnod(i,1),1)=Vnod(i,2);

end

Solution

Vo=A\B; 

disp(Vo);                                                

Getting Started with Simulink: Simple Power calculations

Getting Started with Simulink: Simple Power calculations

In this small example, we will calculate the instantaneous power along with average rms power and voltage using simulink when equation of V and I are given.

Problem:
if v(t)=10cos(120πt+30) and i(t)=6cos(120πt+60).Determine the average power and rms value of v(t).

Formulas Used:
As we know instantaneous power is given as
\[P(t)=v(t) i(t) \]
and average power is given as
\[P_{avg}=\frac{1}{T} \int_0^T{P(t)dt} \]
similarly rms voltage is given as
\[V_{rms}=\sqrt {\frac{1}{T} \int_0^T{V^2(t)dt}} \]
where T is an integer multiple of timeperiod.

Source Creation:
First we will create two sin sources V and I and adjust their parameter to match the problem
V Source:
 Amplitude: 10
  freq:          60
  phase:       pi/6
 sample time: 0.01

I Source:
 Amplitude:  6
  freq:          60
  phase:       pi/3
 sample time: 0.01

RMS Voltage
we will simple follow the equation of Vrms for the calculation.
1. first put a math function block, double click on it and change the equation as u^2
    connect this block's input from the output of V source.
    Output of this block represent V^2
 2. Now put a integrator block and connect the input of it from V^2. The output represents \[ \int_0^T{V^2(t)dt} \]
 3. Now put a divider block and connect the output of integrator to its input.
     Output of this block represents                    \[\frac{1}{T} \int_0^T{V^2(t)dt} \]
 4. Put a math function block, double click on it and change the equation as sqrt. Now connect this to
     the output of step3.
 5. Put a display block and connect the step4 output to it.

Fig1: Simulink model for RMS Voltage calculation

Instantaneous Power
Calculation of Inst, Power is very simple. Just multiply the V and I source and put a scope block to display it.

Fig2: Inst. Power

 Average Power
 Again Average Power calculation is simple and follows from RMS Voltage calculation. Put an integrator and divide block after the instant. power and connect it with the display block.
.

Multipages GUI forms: combining from muliple GUI m files : Links approach

We have learned to create a multipage form single GUI. This is not a very clean approach, because all the panels are on same position and it is sometime very hard to edit one panel later on. This example shows you to create different GUI for different  pages/slides and create a master gui to control them.

(See the previous Note on Multipage GUIs, at http://www.facebook.com/notes/matlab-by-examples-book/gui-working-on-multislides-form/155932614460266)

Handles

 Remember:
  Each GUI has a structure called handles. To see or edit this  structure , you can call guidata.

GET HANDLES

h=guidata(gui_reference);

SET HANDLES

guidata(gui_reference,h);

gui_reference is a double number which works as a pointer to the gui. To store this pointer as a variable, you need to call gui with an output

gui_reference=gui1;

Note that gui's handles.output also contain this pointer. You can use this output when you are writing code inside that gui.

Creating Slide Pages
 First of All create two guis both (GUI1.m and GUI2.m) with two edit box edit1,edit2 and two pushbutton pushbutton1 with title Next and second  pushbutton pushbutton2 with title Previous.

Now for navigation we will add too fields to each gui, next and prev.
So first in OpeningFcn add this two lines both in gui1 and gui2

function gui1_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject    handle to figure

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

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

% varargin   command line arguments to a1 (see VARARGIN)

% Choose default command line output for a1

handles.output = hObject;
handles.next = 1 ;
handles.prev = 1 ;

% Update handles structure

guidata(hObject, handles);

Now go to pushbutton1_callback and add these line for displaying the next slide and disappearing itself

set(handles.next,'Visible','On');
set(handles.output,'Visible','Off');

for pushbutton2_callback do this to display previous slide.

set(handles.prev,'Visible','On');
set(handles.output,'Visible','Off');

Creating Master Slide
Now create a blank GUI with a pushbutton titled "Start"  and save it as gui0.fig. This will automatically create a gui0.m file.

Now in this file, go to guiOpeningFcn
and you will find the followin code there

function gui0_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject    handle to figure

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

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

% varargin   command line arguments to gui0 (see VARARGIN)

% Choose default command line output for gui0

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

Now edit this code to call two slides and save their references to it.

% Choose default command line output for a1

handles.output = hObject;


handles.s1= gui1;
handles.s2= gui2;


h1=guidata(handles.s1);
h1.next = handles.s2;
h1.prev = hObject;
guidata(handles.s1,h1);


h2=guidata(handles.s2);
h2.prev = handles.s1;
h2.next = hObject;
guidata(handles.s2,h2);

% Update handles structure

guidata(hObject, handles);
handles.output

%set(handles.output,'Visible','off');

set(handles.s1,'Visible','off');
set(handles.s2,'Visible','off');guidata(hObject, handles);

As you can see, I have called gui1 and gui2 and saved there reference or pointers to handles.s1 and handles.s2. Now for navigation, we added the s2 to the next of s1(slide 1) and mastergui as previous of s1. Similarly for s2.

Now we need to code the start button to hide the master page and start the silde1.
Here it is .. simple enough now..

set(handles.output,'Visible','off');
set(handles.s1,'Visible','on');
set(handles.s2,'Visible','off');

remember that handles.output indicate the pointer to gui itself.


Final Touch: Processing your data
Now once you are done , you return to master page again because next slide of the s2 is masterpage.
Here you need to collect all the data and do whatever you want. As an example I will add all the 4 number together which were written in 4 editboxes. For easeness, Put one more button titled "Submit" and one Editbox titled Result on masterpage and write the calculation code in Submit_callback. Needless to say that you need to press this submit button once you are back on masterpage to do the final editing.

Here is the calculation code for your reference

h1=guidata(handles.s1);
a1=get(h1.'edit1','String');
a2=get(h1.'edit2','String');
h2=guidata(handles.s2);
a3=get(h2.'edit1','String');
a4=get(h2.'edit2','String');
a=str2num(a1)+str2num(a2)+str2num(a3)+str2num(a4);
set(handles.result,'String',a);

Here you are done. The benefit of the master page is that you can rearrange your slides very easily and each gui is independent.

Generation of Surface by Radar Data: 1 D Example

This Example shows how you can generate a surface from radar data in 1D plane.

Setup
In the setup, we take a 1D example where we rotate the radar direction in a plane from theta=0 to 180 and record the time t taken by signal to come back. Lets assume the velocity is v, then the distance r is
 \[r=v \times t/2\]

So in the setup if we take steps of deltatheta in rotation , we have thetamat matrix denoting different rotation
position and correspondingly we have tmat matrix from which calculate the rmat matrix.

Fig 1 Surface generation from Radar

Now calculating the surface(line here) is pretty easy.

x=rcos(theta)
y=-rsin(theta)

assuming radar is at (0,0)

Generation of Radar Data
Since we are not doing experiment here we need to generate radar data by simulation. we will generate random data using rand function.

thetamat=0:0.01:2*pi;
%generation of Radar Data(use the real recorded data here
tmat=rand(size(thetamat));
%for more smooth surface remove the comments below
%windowSize = 5;
%tmat=filter(ones(1,windowSize)/windowSize,1,1+tmat*0.3);

Generation of Surface
Here is the code:

v=1000;
rmat=v*tmat/2;

xmat=rmat.*cos(thetamat);
ymat=-rmat*.sin(thetamat);

plot(xmat,ymat);

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