Ploting Axes in Matlab

1. Initial

MATLAB has extensive facilities for displaying vectors and matrices as graphs, as well as
annotating and printing these graphs. This section describes a few of the most important graphics functions and provides examples of some typical applications.

The plot function has different forms, depending on the input arguments. If y is a vector, plot(y) produces a piecewise linear graph of the elements of y versus the index of the elements of y. If you specify two vectors as arguments, plot(x,y) produces a graph of y versus x.

For example, to plot the value of the sine function from zero to 2, use

t = 0:pi/100:2*pi;
y = sin(t);
plot(t,y)

Multiple x-y pairs create multiple graphs with a single call to plot. MATLAB automatically cycles through a predefined (but user settable) list of colors to allow discrimination between each set of data. For example, these statements plot three related functions of t, each curve in a separate distinguishing color:

y2 = sin(t-.25); y3 = sin(t-.5); plot(t,y,t,y2,t,y3)

It is possible to specify color, linestyle, and markers, such as plus signs or circles, with:

plot(x,y,'color_style_marker')

color_style_marker is a 1-, 2-, or 3-character string (delineated by single quotation marks) constructed from a color, a linestyle, and a marker type:

• Color strings are 'c', 'm', 'y', 'r', 'g', 'b', 'w', and 'k'. These correspond to cyan, magenta, yellow, red, green, blue, white, and black.
• Linestyle strings are '-' for solid, '- -' for dashed, ':' for dotted, '-.' for dash-dot, and 'none' for no line.
• The most common marker types include '+', 'o', '*', and 'x'.

For example, the statement:

plot(x,y,'y:+')

plots a yellow dotted line and places plus sign markers at each data point. If you specify a marker type but not a linestyle, MATLAB draws only the marker.

2. Figure Windows

The plot function automatically opens a new figure window if there are no figure windows already on the screen. If a figure window exists, plot uses that window by default. To open a new figure window and make it the current figure, type

figure

To make an existing figure window the current figure, type

figure(n)

where n is the number in the figure title bar. The results of subsequent graphics commands are displayed in this window.


3. Adding Plots to an Existing Graph

The hold command allows you to add plots to an existing graph. When you type

hold on

MATLAB does not remove the existing graph; it adds the new data to the current graph, rescaling if necessary. For example, these statements first create a contour plot of the peaks function, then superimpose a pseudocolor plot of the same function:

[x,y,z] = peaks; contour(x,y,z,20,'k') hold on pcolor(x,y,z) shading interp

The hold on command causes the pcolor plot to be combined with the contour plot in one figure.

4. Subplots

The subplot function allows you to display multiple plots in the same window or print them on the same piece of paper. Typing

subplot(m,n,p)

breaks the figure window into an m-by-n matrix of small subplots and selects the pth subplot for the current plot. The plots are numbered along first the top row of the figure window, then the second row, and so on. For example, to plot data in four different subregions of the figure window,

t = 0:pi/10:2*pi;
[X,Y,Z] = cylinder(4*cos(t));
subplot(2,2,1); mesh(X)
subplot(2,2,2); mesh(Y)
subplot(2,2,3); mesh(Z)
subplot(2,2,4); mesh(X,Y,Z)

5. Imaginary and Complex Data

When the arguments to plot are complex, the imaginary part is ignored except when plot is given a single complex argument. For this special case, the command is a shortcut for a plot of the real part versus the imaginary part. Therefore,

plot(Z)

where Z is a complex vector or matrix, is equivalent to

plot(real(Z),imag(Z))

For example

t = 0:pi/10:2*pi;
plot(exp(i*t),'-o')

draws a 20-sided polygon with little circles at the vertices.

6. Controlling Axes

The axis function has a number of options for customizing the scaling, orientation, and aspect ratio of plots.

Ordinarily, MATLAB finds the maxima and minima of the data and chooses an appropriate plot box and axes labeling. The axis function overrides the default by setting custom axis limits,

axis([xmin xmax ymin ymax])

axis also accepts a number of keywords for axes control.   For example

axis square

makes the entire x-axes and y-axes the same length and

axis equal

makes the individual tick mark increments on the x- and y-axes the same length. So

plot(exp(i*t))

followed by either axis square or axis equal turns the oval into a proper circle.

axis auto

returns the axis scaling to its default, automatic mode.

axis on

turns on axis labeling and tick marks.

axis off

turns off axis labeling and tick marks.

The statement

grid off

turns the grid lines off and

grid on

turns them back on again.

7. Axis Labels and Titles

The xlabel, ylabel, and zlabel functions add x-, y-, and z-axis labels. The title function adds a title at the top of the figure and the text function inserts text anywhere in the figure. A subset of Tex notation produces Greek letters, mathematical symbols, and alternate fonts. The following example uses \leq for , \pi for , and \it for italic font.

t = -pi:pi/100:pi;
y = sin(t);
plot(t,y)
axis([-pi pi -1 1])
xlabel('-\pi \leq {\itt} \leq \pi')
ylabel('sin(t)')
title('Graph of the sine function')
text(1,-1/3,'\it{Note the odd symmetry.}')

8. Mesh and Surface Plots

MATLAB defines a surface by the z-coordinates of points above a grid in the x-y plane, using straight lines to connect adjacent points. The functions mesh and surf display surfaces in three dimensions. mesh produces wireframe surfaces that color only the lines connecting the defining points. surf displays both the connecting lines and the faces of the surface in color.

9. Visualizing Functions of Two Variables

To display a function of two variables, z = f (x,y), generate X and Y matrices consisting of repeated rows and columns, respectively, over the domain of the function. Then use these matrices to evaluate and graph the function. The meshgrid function transforms the domain specified by a single vector or two vectors x and y into matrices X and Y for use in evaluating functions of two variables. The rows of X are copies of the vector x and the columns of Y are copies of the vector y.

To evaluate the two-dimensional sinc function, sin(r)/r, between x and y directions:

[X,Y] = meshgrid(-8:.5:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R)./R;
mesh(X,Y,Z)




In this example, R is the distance from origin, which is at the center of the matrix. Adding eps avoids the indeterminate 0/0 at the origin.

10. Images

Two-dimensional arrays can be displayed as images, where the array elements determine brightness or color of the images. For example,

load durer
whos

shows that file durer.mat in the demo directory contains a 648-by-509 matrix, X, and a 128-by-3 matrix, map. The elements of X are integers between 1 and 128, which serve as indices into the color map, map. Then

image(X)
colormap(map)
axis image

reproduces Dürer's etching shown at the beginning of this book. A high resolution scan of the magic square in the upper right corner is available in another file. Type

load detail

and then use the uparrow key on your keyboard to reexecute the image, colormap, and axis commands. The statement

colormap(hot) 

adds some twentieth century colorization to the sixteenth century etching.


11. Printing Graphics

The Print option on the File menu and the print command both print MATLAB figures. The Print menu brings up a dialog box that lets you select common standard printing options. The print command provides more flexibility in the type of output and allows you to control printing from M-files. The result can be sent directly to your default printer or stored in a specified file. A wide variety of output formats, including PostScript, is available.

For example, this statement saves the contents of the current figure window as color Encapsulated Level 2 PostScript in the file called magicsquare.eps:

print -depsc2 magicsquare.eps

It's important to know the capabilities of your printer before using the print command. For example, Level 2 Postscript files are generally smaller
and render more quickly when printing than Level 1 Postscript. However, not all PostScript printers support Level 2, so you need to know what
your output device can handle. MATLAB produces graduated output for surfaces and patches, even for black and white output devices. However,
lines and text are printed in black or white.

Simulink Introduction

Simulink (Simulation and Link) is an extension of MATLAB. It works with MATLAB to offer modeling, simulating, and analyzing of dynamical systems under a graphical user interface (GUI) environment. The construction of a model is simplified with click-and-drag mouse operations. Simulink includes a comprehensive block library of toolboxes for both linear and nonlinear analysis. Models are hierarchical, which allow using both top-down and bottom-up approaches. As Simulink is an integral part of MATLAB, it is easy to switch back and forth during the analysis process and thus, the user may take full advantage of features offered in both environments. This tutorial presents the basic features of Simulink and is focused on control systems. Very useful for engineering field. This tutorial has been written for Simulink v.5 and v.6 It is upward and backward version compatible.

Getting Started

To start a Simulink session, you'd need to bring up Matlab program first.

From Matlab command window, enter:

>> simulink

Alternately, you may click on the Simulink icon located on the toolbar as shown:

Simulink's library browser window like one shown below will pop up presenting the block set for model construction.
To see the content of the blockset, click on the "+" sign at the beginning of each toolbox. To start a model click on the NEW FILE ICON as shown in the screenshot above. Alternately, you may use keystrokes CTRL+N. A new window will appear on the screen. You will be constructing your model in this window. Also in this window the constructed model is simulated. A screenshot of a typical working (model) window is shown below:


To become familiarized with the structure and the environment of Simulink, you are encouraged to explore the toolboxes and scan their contents. You may not know what they are all about at first, but perhaps you could catch on the organisation of these toolboxes according to their categories. For instance, you may see that the Control System toolbox consists of the Linear Time Invariant (LTI) system library and the MATLAB functions can be found under Function and Tables of the Simulink main toolbox. A good way to learn Simulink (or any computer program in general) is to practice and explore. Making mistakes is part of the learning curve. So, fear not you should be!
A simple model is used here to introduce some basic features of Simulink. Please follow the steps below to construct a simple model.

STEP 1: CREATING BLOCKS.

From BLOCK SET CATEGORIES section of the SIMULINK LIBRARY BROWSER window, click on the "+" sign next to the Simulink group to expand the tree and select (click on) Sources.

















A set of blocks will appear in the BLOCKSET group. Click on the Sine Wave block and drag it to the workspace window (also known as model window).


Now you have established a source of your model.

NOTE: It is advisable that you save your model at some point early on so that if your PC crashes you wouldn't loose too much time reconstructing your model.

I am going to save this model under the filename: "sinexample1". To save a model, you may click on the floppy diskette icon . or from FILE menu, select Save or using keystrokes CTRL+S. All Simulink model file will have an extension ".mdl". Simulink recognises file with .mdl extension as a simulation model (similar to how MATLAB recognises files with the extension .m as an MFile).
Continue to build your model by adding more components (or blocks) to your model window. We'll continue to add a Scope from Sinks library, an Integrator block from Continuous library, and a Mux block from Signal Routing library.

NOTE: If you wish to locate a block knowing its name, you may enter the name in the SEARCH WINDOW (at Find prompt) and Simulink will bring up the specified block.

To move the blocks around, simply click on it and drag it to a desired location. Once you've dragged over all necessary blocks, the workspace window should consist of the following components:


You may remove (delete) a block by simply clicking on it once to turn on the "select mode" (with four corner boxes) and use the DEL key or keys combination CTRL-X.

STEP 2: MAKING CONNECTIONS

To establish connections between the blocks, move the cursor to the output port represented by ">" sign on the block. Once placed at a port, the cursor will turn into a cross "+" enabling you to make connection between blocks.
A sine signal is generated by the Sine Wave block (a source) and is displayed by the scope. The integrated sine signal is sent to scope for display along with the original signal from the source via the Mux, whose function is to mutiplex signals in form of scalar, vector, or matrix into a bus.

STEP 3: RUNNING SIMULATION

You now may run the simulation of the simple system above by clicking on the play button . Alternately, you may use keystrokes CTRL+T, or choose Start submenu (under Simulation menu).
Double click on the Scope block to display of the scope.


To view/edit the parameters, simply double click on the block of interest.

Handling of Blocks and Lines

The table below describes the actions and the corresponding keystrokes or mouse operations (Windows versions).

Actions	Keystokes or Mouse Actions
Copying a block from a library	Drag the block to the model window with the left button on the mouse OR use select the COPY and PASTE from EDIT menu.
Duplicating blocks in a model	Hold down the CTRL key and select the block with the left mouse drag the block to a new location.
Display block's parameters	Double click on the block
Flip a block	CTRL-F
Rotate a block (clockwise 90 deg @ each keystroke)	CTRL-R
Changing blocks' names	Click on block's label and position the cursor to desired place.
Disconnecting a block	hold down the SHIFT key and drag the block to a new location
Drawing a diagonal line	hold down the SHIFT key while dragging the mouse with the left button
Dividing a line	move the cursor to the line to where you want to create the vertex and use the left button on the mouse to drag the line while holding down the SHIFT key

Annotations

To add an annotation to your model, place the cursor at an unoccupied area in your model window and double click (left button). A small rectangular area will appear with a cursor prompting for your input.

To delete an annotation, hold down the SHIFT key while selecting the annotation, then press the DELETE or BACKSPACE key. You may also change font type and colour from FORMAT menu. Make sure that the block is selected before making the change.

Hope this useful. Thanks for visiting.

Introduction to Matlab

Matlab is both a powerful computational environment and a programming language that easily handles matrix and complex arithmetic. It is a large software package that has many advanced features built-in, and it has become a standard tool for many working in science or engineering disciplines. Among other things, it allows easy plotting in both two and three dimensions.

Matlab has two different methods for executing commands: interactive mode and batch mode. In interactive mode, commands are typed (or cut-and-pasted) into the 'command window'. In batch mode, a series of commands are saved in a text file (either using Matlab's built-in editor, or another text editor such as Emacs) with a '.m' extension. The batch commands in a file are then executed by typing the name of the file at the Matlab command prompt. The advantage to using a '.m' file is that you can make small changes to your code (even in different Matlab sessions) without having to remember and retype the entire set of commands. Also, when using Matlab's built-in editor, there are simple debugging tools that can come in handy when your programs start getting large and complicated. More on writing .m files later.

Scalar Variables and Arithmetic Operators

Scalar variables are assigned in the obvious way:

>> x = 7
x =
7

The >> is Matlab's command prompt. Notice that Matlab echos back to you the value of the assignment you have made. This can be helpful occasionally, but will generally be unwieldy when working with vectors. By placing a semicolon at the end of a statement, you can suppress this verbose behavior.

>> x = 7;

Variables in Matlab follow the usual naming conventions. Any combination of letters, numbers and underscore symbols ('_') can be used, as long as the first character is a letter. Note also that variable names are case sensitive ('X' is different from 'x').

All of the expected scalar arithmetic operators are available:

>> 2*x
ans =
14

>> x^2
ans =
49

Notice that the multiply operation is not implied as it is in some other computational environments and the '*' operator has to be specified (typing '>>2x' will result in an error). This is also a good time to point out that Matlab remembers the commands that you are entering and you can use the up-arrow button to scroll through them and edit them for re-entry. This is very handy, especially when you are repeatedly making minor changes to a long command line.

There are a few predefined variables in Matlab that you will use often: pi, i, and j. Both i and j are defined to be the square-root of -1 and pi is defined as the usual 3.1416... . These variables can be assigned other values, so the statement

>> pi = 4;

is valid. Care is needed to avoid changing a predefined variable and then forgetting about that change later. That may seem unreasonable with pi, but i and j are both natural variables to use as indices and they are often changed. A useful command is the 'clear', or 'clear all' command. Typing either of these at the command prompt will remove all current variables from Matlab's memory and reset predefined variables to their original values. The clear command can also remove one specific variable from memory. The command

>> clear x

would only remove the variable named 'x' from Matlab's memory.

Matrix Variables and Arithmetic Operators

Another useful Matlab command is 'whos', which will report the names and dimensions of all variables currently in Matlab's memory. After typing the command at the prompt we see:

>> whos
Name Size Bytes Class

x 1x1 8 double array

Grand total is 1 elements using 8 bytes

The important thing to note right now is that the size is given as '1x1'. Matlab is really designed to work with vectors and matrices, and a scalar variable is just a special case of a 1x1 dimensional vector. To assign a vector containing the first 5 integers to the variable x, we could type this command:

>> x = [1 2 3 2^2 2*3-1]
x =
1 2 3 4 5

We won't have much occasion to operate on matrices that are higher dimension, but if you wanted to create a 2-D matrix you could use a command something like:

>> A = [1 2 3; 4 5 6]
A =
1 2 3
4 5 6

To create larger vectors than the toy examples above (say, the integers up to 100), we would need to type a lot of numbers. Not surprisingly, there are easier ways built in to create vectors. To create the same 5-element vector we did above, we could also type:

>> x = [1:5]
x =
1 2 3 4 5

The colon operator creates vectors of equally spaced elements given a beginning point, and maximum ending point and the step size in between elements. Specifically, [b:s:e] creates the vector [b b+s b+2*s ... e]. If no step size is specified (as in the example above), a step of 1 is assumed. So, the command [1:2:10] would create the vector of odd integers less than 10, [1 3 5 7 9] and the command [1:3:10] would create the vector of elements [1 4 7 10].

Let's create the vector of odd elements mentioned above:

>> x_odd = [1:2:10]
x_odd =
1 3 5 7 9

You can access any element by indexing the vector name with parenthesis (indexing starts from one, not from zero as in many other programming languages). For instance, to access the 3rd element, we would type:

>> x_odd(3)
ans =
5

We can also access a range of elements (a subset of the original vector) by using the colon operator and giving a starting and ending index.

>> x_odd(2:4)
ans =
3 5 7

If we want to do simple arithmetic on a vector and a scalar, the expected things happen,

>> 3+[1 2 3]
ans =
4 5 6

>> 3*[1 2 3]
ans =
3 6 9

and the addition or subtraction of matrices is possible as long as they are the same size:

>> [1 2 3]+[4 5 6]
ans =
5 7 9

The operators '*' and '/' actually represent matrix multiplication and division which is not typically what we will need in this course. However, a common task will be to form a new vector by multiplying (or dividing) the elements of two vectors together, and the special operators '.*' and './' serve that purpose.

>> [1 2 3].*[4 5 6]
ans =
4 10 18

>> [1 2 3]./[4 5 6]
ans =
0.2500 0.4000 0.5000

Beware that the operator '^' is a shortcut for repeated matrix multiplications ('*'), whereas the operator '.^' is the shortcut for repeated element-by-element multiplications ('.*'). So, to square all of the elements in a vector, we would use

>> [1 2 3].^2 ans = 1 4 9

Built-in Commands

Matlab has many built-in commands to do both elementary mathematical operations and also complex scientific calculations. The 'help' command will be useful in learning about built-in commands. By typing help at the command prompt, you are given a list of the different categories that Matlab commands fall into (e.g., general, elementary matrix operations, elementary math functions, graphics, etc.). Notice that there are probably even specific toolboxes included with your Matlab package for performing computations for disciplines such as signal processing. To see all of the commands under a certain topic, type 'help topic'. To get a description of a specific command, type 'help command'.

We'll look at a couple of example commands. First, the square-root command, sqrt. To calculate the square root of a number, type:

>> sqrt(2)
ans =
1.4142

Again, to illustrate that Matlab understands complex numbers, we can have it calculate the root of a negative number:

>> sqrt(-9)
ans =
0 + 3.0000i

Many Matlab commands that operate on scalars will also operate element-by-element if you give it a vector. For instance:

>> sqrt([1 2 4 6])
ans =
1.0000 1.4142 2.0000 2.4495

>> sqrt([1:8])
ans =
1.0000 1.4142 1.7321 2.0000 2.2361 2.4495 2.6458 2.8284

Another useful command is sin function, which also operates on a vector.

>> sin([pi/4 pi/2 pi])
ans =
0.7071 1.0000 0.0000

It is important to realize that digital signals are really just a collection of points and can be represented by a vector. Because Matlab allows functions like sin and sqrt (as well as many others) to operate on vectors, many of the signal definitions and calculation we would like to perform are very straight-forward. We will illustrate this a little more explicitly in the next section.

Signals, Plotting and Batch Operation

We will now use the elementary tools in Matlab to build a basic signal we have talked about in class. All of the commands given in this section are already written in the file sinplot.m. They can be executed in batch by typing sinplot at the command prompt. If the file is not found, use the 'current directory' browser at the top of the Matlab window to make the current directory the one where the file lives.

First, we will define a signal which is a 2 Hz sinewave over the interval [0,1] seconds:

>> t = [0:.01:1]; % independent (time) variable
>> A = 8; % amplitude

>> f_1 = 2; % create a 2 Hz sine wave lasting 1 sec
>> s_1 = A*sin(2*pi*f_1*t);

Anything following a '%' symbol is a comment and will be ignored by Matlab. A 4 Hz sinewave with the same amplitude will also be defined:

>> f_2 = 4; % create a 4 Hz sine wave lasting 1 sec
>> s_2 = A*sin(2*pi*f_2*t);

There are a few commands necessary to do basic plotting of these signals. The figure command will bring up a new figure window and the close command closes a specific window ('close all' closes all open windows). The subplot(r, c, p) command allows many plots to be 'tiled' into r rows and c columns on one figure window. After this command is issued, any succeeding commands will affect the p^th tile. The plot(x,y) command will then plot the equal-length vectors x and y on the appropriate axis and in the obvious way. The two sinewaves created above are then plotted along with the sum of the two signals in separate tiles with the commands:

%plot the 2 Hz sine wave in the top panel
figure
subplot(3,1,1)
plot(t, s_1)
title('2 Hz sine wave')
ylabel('Amplitude')

%plot the 4 Hz sine wave in the middle panel
subplot(3,1,2)
plot(t, s_2)
title('4 Hz sine wave')
ylabel('Amplitude')

%plot the summed sine waves in the bottom panel
subplot(3,1,3)
plot(t, s_1+s_2)
title('Summed sine waves')
ylabel('Amplitude')
xlabel('Time (s)')

The commands title, xlabel, and ylabel all take strings as arguments and apply the appropriate label to the top, x-axis or y-axis of the current subplot, respectively.

Normally, whenever you have plotted a function on a subplot and you try to plot another function, Matlab will erase the subplot and start over. Often you will want to plot two signals on top of each other for comparison. This can be done by using the 'hold', which toggles the hold property of a subplot on or off. In this next example, the 2 original sinewaves are created again, this time using the built-in exponential function exp, the built-in function imag and Euler's relation. The imag function just returns the imaginary part of a complex number. These statements accomplish the same thing that the sin function did earlier, and are only used to illustrate again that Matlab handles all of the complex arithmetic we need very naturally. In this example, the two sinewaves are plotted on top of each other, and the sum is plotted again in its own subplot. Notice that the plot command can also take a third argument that specifies the color and linestyle of the resulting plot. There are many options here, and they can be viewed by typing help plot at the prompt.

% create the same sine waves using Euler's relations
s_3 = (exp(j*2*pi*f_1*t) - exp(-j*2*pi*f_1*t))/(2*j);
s_4 = imag(exp(j*2*pi*f_2*t));

% plot the 2 and 4 Hz waves together in the same panel
figure
subplot(2,1,1)
plot(t,s_3, 'b-')
hold on
plot(t, s_4, 'r--')
ylabel('Amplitude')

% again plot the sum in the bottom panel
subplot(2,1,2)
plot(t, s_3+s_4)
ylabel('Amplitude')
xlabel('Time (s)')

Basic Programming (For Loops)

Though we haven't mentioned much of this, Matlab is also actually a programming language with all of the usual iterative and conditional execution abilities. Because it is matrix based, many operations that you would normally use a loop for (say, if you were programming in C) can be done much faster using Matlab's built-in matrix operations. But, loops are sometimes unavoidable when you have repeated operations to perform on a list of parameters, and we will illustrate a simple one with an example here. All of the commands given in this section can be found in the file sumsinplot.m, and can be executed in batch by typing sumsineplot at the command prompt.

By far, the most common looping apparatus (in Matlab, at least) is the for loop. The basic for loop has the following structure:

for variable = values
statements
end

where variable in this case is generally some sort of counter or indexing variable. This is probably best explained through an illustration. In this example, we will generate sinwaves with odd integers for frequencies (1, 3, 5, 7, 9, 11) and cumulatively sum them together. After each new sinewave is added to the others, the result is plotted in its own subplot. To be explicit, a 1 Hz sinewave is generated and plotted. Then a 3 Hz sinewave is generated, added to the 1 Hz sinewave and plotted. A 5 Hz sinewave is generated, added to the sum of the 1st two signals and plotted, etc. The only new things in this example are the use of the for loop, and the use of the zeros command. The zeros(n,m) command returns a vector of dimension (n x m), initialized to all zeros.

%plot cumulatively summed sine waves with odd integer frequencies scaled by 1/f

>> t = 0:.001:1; % independent (time) variable
>> A = 1; % amplitude
>> f = 1:2:11; % odd frequencies
>> s = zeros(1,101); % empty 'signal' vector to start with

>> figure

% for each frequency, create the signal and add it into s(t)
>> for count=1:6
>> s = s + A*sin(2*pi*f(count)*t)/f(count);
>> subplot(6,1,count)
>> plot(t,s);
>> ylabel('Amplitude')
>> end

>> subplot(6,1,1)
>> title('Cumulatively summing sine waves')
>> subplot(6,1,6)
>> xlabel('Time (s)')

Imfilter function for Nonlinear Operation

Some nonlinear image processing operations can be expressed in terms of linear filtering. When this is true, it often provides a recipe for a speedy MATLAB implementation. Today I'll show two examples: Computing the local standard deviation, and computing the local geometric mean.

The local standard deviation operator is often used as a measure of "busy-ness" throughout the image. For each pixel, the standard deviation of that pixel's neighbors is computed:

sqrt(sum((x - xbar).^2));

Mathematically, the summation can be rewritten as:

sum(x.^2) - sum(x).^2/N

Both of these local sums can be computed by using imfilter with an all-ones filter.

I = im2double(imread('cameraman.tif')); imshow(I);

To compute the sum(x.^2) term, we square (elementwise) the input to imfilter.

h = ones(5,5);
term1 = imfilter(I.^2, h);
imshow(term1, [])

Notice the dark band around the edge. This is because imfilter zero-pads by default. We might want to use the 'symmetric' option instead.

term1 = imfilter(I.^2, h, 'symmetric');
imshow(term1, [])

To compute the sum(x).^2 term, we square the output of imfilter.

term2 = imfilter(I, h, 'symmetric').^2 / numel(h);
imshow(term2, [])

Then we subtract the second term from the first and take the square root.

local_std = sqrt(term1 - term2);  % scale factor omitted
imshow(local_std, [])

Cautionary Notes

1. The procedure shown here is not always a numerically sound way to compute the standard deviation, because it can suffer from both overflow (squared terms) and underflow (cancellation in the subtraction of large numbers). For typical image pixel values and window sizes, though, it works reasonably well.

2. Round-off error in the computation of term1 and term2 can sometimes make (term1 - term2) go slightly negative, resulting in complex outputs from square root operator. Avoid this problem by using this code:

local_std = sqrt(max(term1 - term2, 0));

Recent releases of the Image Processing Toolbox include the function stdfilt, which does all this work for you.
The local geometric mean filter multiplies together all the pixel values in the neighborhood and then takes the N-th root, where N is the number of pixels in the neighborhood. The geometric mean filter is said to be slightly better than the arithmetic mean at preserving image detail.
Use the old logarithm trick to express the geometric mean in terms of a summation. Then imfilter can be used to compute the neighborhood summation, like this.

geo_mean = imfilter(log(I), h, 'replicate');
geo_mean = exp(geo_mean);
geo_mean = geo_mean .^ (1/numel(h));

imshow(geo_mean, [])

This is the complete source code. You can directly copy and paste to an M-File and run it.

%% Nonlinear filtering using imfilter
% Some nonlinear image processing operations can be
% expressed in terms of linear filtering. 
% When this is true, it often provides a recipe
% for a speedy MATLAB implementation. Today I'll
% show two examples:
% Computing the local standard deviation, and
% computing the local geometric mean.

%%
% The local standard deviation operator
% is often used as a measure of "busy-ness"
% throughout the image.  For each pixel,
% the standard deviation of that
% pixel's neighbors is computed:
%
% sqrt( sum ( (x - xbar).^2 ) )
%
% (Scale factors omitted.)
%
% Mathematically, the summation can be
% rewritten as:
%
% sum(x.^2) - sum(x).^2/N
%
% Both of these local sums can be computed
% by using with an all-ones filter.

I = im2double(imread('cameraman.tif'));
imshow(I)

%%
% _Original image credit:
% Massachusetts Institute of Technology_
%
% To compute the sum(x.^2) term, we square
% (elementwise) the input to |imfilter|.

h = ones(5,5);
term1 = imfilter(I.^2, h);
imshow(term1, [])

%%
% Notice the dark band around the edge. 
% This is because |imfilter| zero-pads
% by default.  We might want to use the
% 'symmetric'option instead.

term1 = imfilter(I.^2, h, 'symmetric');
imshow(term1, [])

%%
% To compute the sum(x).^2 term, we square
% the output of |imfilter|.

term2 = imfilter(I,h,'symmetric').^2/numel(h);
imshow(term2, [])

%%
% Then we subtract the second term from
% the first and take the square root.

local_std = sqrt(term1-term2); % scale factor omitted
imshow(local_std, [])

%%
% _Cautionary notes_
%
% * The procedure shown here is not always a numerically
% sound way to compute the standard deviation,
% because it can suffer from both overflow
% (squared terms) and underflow (cancellation in the
% subtraction of large numbers). For typical image pixel
% values and window sizes, though, it works reasonably
% well.
%
% * Round-off error in the computation of term1 and
% term2 can sometimes make (term1 - term2) go slightly
% negative, resulting in complex outputs from square root
% operator. Avoid this problem by using this code:

local_std = sqrt(max(term1 - term2, 0));

%%
% Recent releases of the Image Processing Toolbox
% include the function which does all this work for you.

%%
% The local geometric mean filter multiplies together
% all the pixel values in the neighborhood
% and then takes the N-th root, where N is the number
% of pixels in the neighborhood. The geometric mean
% filter is said to be slightly better than
% the arithmetic mean at preserving image detail.
%
% Use the old logarithm trick to express the geometric
% mean in terms of a summation.  Then |imfilter| can be
% used to compute the neighborhood summation, like this.

geo_mean = imfilter(log(I), h, 'replicate');
geo_mean = exp(geo_mean);
geo_mean = geo_mean .^ (1/numel(h));

imshow(geo_mean, [])

Copyright 2008 Mathworks, Inc