Tuesday 22 February 2011

MATHWORkS MATLAB TUTORIAL AND EXERCISE


MATHWORKS MATLAB: 

Exercise #1: Non-oriented, Center-Surround Receptive Field Model.

Goal: To implement and test a simple RF model in Matlab.

1. MATHWORKS MATLAB Basics

When you start Matlab, it should open with a Command Window.  This window allows you to enter individual commands at the prompt (>>).  Commands can define variables.  For example:

>> x = zeros(100);

creates a 2D array (100x100) and fills it with zeros.  Variables are initialized whenever they appear on the left hand side of an expression.  You do not need to specify the storage type or size, as in other languages.  Variables must be initialized before they are used on the right hand side of an expression.  When variables are initialized they get stored in the Workspace.  You can view the variables in the Workspace by typing:

>> who

You can see the contents of a single variable by typing the variable name at the prompt without a semicolon at the end:

>> x

Commands can also call Matlab functions, e.g.

>> plot(x);

To find out more about the function “plot”, type:

>> help plot
           
If you want to create a program of more than a few lines, it soon becomes very tedious to enter and re-enter commands one at a time.  Therefore, Matlab allows you to create M-files.  M-files are plain text files that contain Matlab expressions and function calls.  M-file names end with the extension “.m”.  Usually, the Command Window has a button that will invoke an M-file editor.  You can type multiple commands into the editor window and save them.  Usually, there will be a button in the editor window that allows you to save and execute the file in one shot.  Otherwise, you  can execute the M-file by typing its name, sans extension, at the command prompt.  E.g. if you created the file “myprog.m”, you would just type:

>> myprog


Important note:  Each time you run an M-file, the variables in the workspace are not cleared.  This can lead to a situation where you spend several hours working on a program and it seems to run fine.  Then you come back to it later, after starting a fresh Matlab session, and you get a bunch of error messages.  What happened?  Well, your program was probably using variables that were in the workspace but not actually initialized by the program itself.  There is a simple way to avoid this.  Start every M-file with the following command:

clear;

BTW, lines in M-files that start with % are treated as comments.

2. A simple RF model.

Retinal ganglion cells and LGN relay cells have, to a first approximation, circularly symmetric receptive fields.  This makes them sensitive to light-dark contrast, but not to contour orientation.  Hence, they are called “non-oriented” or “isotropic” receptive fields.  We will start with a one-dimensional model based on a difference-of-Gaussians (DOG).

Q: Why is one spatial dimension sufficient to model these RFs?

The equation for the RF is the following:

                              

Q: Which terms in this equation correspond to the amplitude, peak position, and spread (width) of the two Gaussians?
                                        
To test the response of this model, we will construct a stimulus that is described by a cosine light intensity function:

Q: Which terms in this equation correspond to the spatial frequency and phase of the cosine?

To get started, we need to define a vector that specifies the range and grain over which we will compute RF(x) and I(x).  So, open an M-file and enter the following commands:

clear;

xmin      = -1.0;
xmax      = 1.0;
dx        = 0.01;
x         = xmin:dx:xmax;

Now, we’ll define the excitatory and inhibitory parts of the RF separately

ke        = 1.0;
xe        = 0.0;
sige      = 0.125;
RFe       = ke * exp(-(x-xe).^2/sige.^2);

Note: The “.^” is not a typo.  Matlab uses the “.” to distinguish simple math functions from matrix functions.  Hence, x.^2 squares each element of x; whereas x^2 multiplies the vector x with itself resulting in a 2D matrix.  The “.” can precede any math function (*, /, …).

Now the inhibitory part:

ki        = -0.25;
xi        = 0.0;
sigi      = 0.5;
RFi       = ki * exp(-(x-xi).^2/sigi.^2);

The full RF is constructed simply by adding the two parts:

RF        = RFe + RFi;

RFe and RFi are vectors of identical length and can therefore be added to produce another vector, RF.  Notice how we were able to do this without a “for” loop.  That’s the beauty of vectorized programming.

Q: What physiological variable might be represented by the vector “RF”?  Firing rate? Membrane potential?  Sensitivity?  Is it “physiologically plausible” for the RF model to use negative numbers?  How might we be able to construct and inhibitory surround without using negative numbers?

At this point, it should be fairly obvious how to construct the stimulus, but I’ll walk you through it anyway.

ws        = 1.0;
phi       = 0.0;
I         = cos(2 * pi * ws * x + phi);

Notice how pi is already defined.  Gotta love those Matlab guys.

To determine the response of our model RF to the stimulus, we need to compute the inner product of the two vectors, RF and I.  This sounds tricky, but it is really very simple:

resp      = dx*sum(I.*RF);

Notice again the use of the “.”  Why is this not needed between “dx” and “sum”?  Why am I multiplying by “dx” in the first place? (Hint: “resp” is the integral of the function resulting from the multiplication of I and RF.  In this formulation, we are implicitly implementing the rectangular rule for integration.  How would you implement the trapezoid or parabola rules?).

At this point, you might want to start plotting some things.  Try this:

figure(1); clf;
subplot(3,1,1);
plot(x,I);
subplot(3,1,2);
plot(x,RF);
subplot(3,1,3);
plot(resp);

To finish the exercise, try the following on your own:

1. Compute the full frequency response by varying the spatial frequency, ws.  A good way to do this is to declare ws as a vector and then use a “for” loop:

ws        = 0.125.*2.^(0:0.5:6);
resp      = zeros(1,length(ws));
for i = 1:length(ws)
     I         = cos(2 * pi * ws(i) * x + phi);
     resp(i)   = dx * sum(I.*RF);
end;
     plot(ws, resp);

Note: if you’re clever, you can accomplish this without a “for” loop by vectorizing the code.

2. Compute the frequency response separately for the excitatory (RFe) and inhibitory (RFi) parts of the RF.  How does the frequency response of the full RF compare to those of its parts?

3. Vary some of the parameters of the RF model and see how this affects the frequency response.  What constraints on the RF parameters are required so that the frequency response doesn’t become negative?

4. Find a way to test the phase response of the RF model.  Design your own stimuli.  The fun never ends!

No comments:

Post a Comment