INTRODUCTION version 04/06/03
Matlab is a wonderful tool for geophysical computation, data analysis and display. Matlab includes about a thousand routines for higher math and data manipulation. The graphics capability is stunning. Student Matlab, which runs under Windows or Linux, contains most of the same routines found in the full professional version, so that students and faculty can work with Student Matlab on their personal computers and use the programs on the professional Matlab. It has taken me a long time to learn what I know about Matlab, and I am learning more all the time. At MTU use Matlab 6 installed on networked Unix machines, and I will describe Matlab in that setting. For running either the full version or Student Matlab on a PC, the concepts are the same, but the methods of starting Matlab and using the built in editor are slightly different.
You should look elsewhere for a comprehensive introduction to Matlab. For a short read, I suggest the book Kermit, Sigmon, A Matlab Primer. Older editions of Sigmon's notes may be downloaded from the World Wide Web, or you may purchase the most recent version. Two excellent longer books which introduce engineering problem solving are:
Chapman, Stephen J., Matlab Programming for Engineers, and
Palm, William J., Introduction to Matlab 6 for Engineers,
These books do not cover the contouring, colorizing and three dimensional data volume manipulations discussed in the notes you are reading now, but they do contain many common computational engineering applications.
The basic data structure in Matlab is the matrix. In programming languages such as Fortran and Basic, this structure is called an array. Matricies can be one, two or three dimensional. The details of working with matrices will be dealt with in more detail in Part II of these notes. For now I will assume that you know what an array or matrix is and that you know how the cells are addressed. Variable types in Matlab can be integer, floating point, single precision, double precision, character, etc. as in most computer languages. Data types will not be otherwise discussed in these notes because it does not affect the material being presented.
Personal note: I am a geophysicist, not a computer scientist. I started learning computer languages with Fortran II in 1964. I learned Basic and more Fortran in graduate school in the 1970s. I have been learning Matlab on my own over the last ten years or so. Please forgive me if I fail to supply basic information that you need. I will try to improvise such explanations if you ask . Please remember: There is no such thing as a stupid question.
To use these workshop notes interactively on your computer: Open a window running Microsoft Word while Matlab is running in another window. You can copy any Matlab example line or short program from the Word document to the Matlab prompt in the window running Matlab. Thanks to Brian Robinson from the University of Lancaster in England for this tip.
Install Matlab now on your PC by following the directions in the booklets in the box.
PART I. FINDING YOUR WAY AROUND MATLAB
This part of the notes is intended for a short workshop or lab session lasting one to four hours.
Start Matlab 6.5 by clicking on the icon on your main Windows screen. Matlab will be ready to use when the computer displays two right brackets (>>). I suggest that you type
demos ,
and take a tour of Matlab examples. The graphics demos are the most relevant to geophysics, including the "slicer" demo. Note that for many of the demos, the short program that generates the example is given in a window. You can copy the program and paste it in your own window and modify it as you like.
All Matlab commands will execute both from .m files and interactively. To execute an .m file you just type its name in the command window, without the .m extension, and press return. To execute a command interactively, you just type it into the command window and press return. In the interactive mode, you generally will loose what you have done, but there is a diary command which will save the commands to a file. For example, entering diary today.txt
will save the commands that follow to the file today.txt.
You may use the interactive mode to experiment with the commands but you should write short programs with an editor program for class work.
Writing Matlab programs with an editor.
You create a matlab program by pulling down the File tab and selecting , New, m file. Matlab programs must have the extension .m. Saving files, etc. is done with the usual Windows pull down tabs. You can save the files and run them immediately in the Matlab window, by typing the name of the file without the .m suffix.
Where do I keep my matlab programs?
When I start up Matlab6, it automatically keeps programs that I write in c:\matlab6p5\work. In fact, it makes sense to create a special directory for each project, for example c:\Matlag\mlnotes and keep your Matlab programs there. You need to set the "path" for matlab so it will know to look in the special directory. Please read the documentation on the path command. I wrote a one line matlab program to do this. The programs is called fixpath.m and it contains this line:
path(path,':\matlab\mlnotes')
If you put the workshop diskette data into c:\matlab6p5\workshop then you should enter:
path(path,‘c:\matlab6p5\workshop‘)
I keep my matlab programs for my matlab notes in the subdirectory . This line sets my path so that Matlab will look in You need to run this program every time you work with programs in that subdirectory. Then Matlab will know to look for your programs there. If this is confusing, talk to a systems administrator or other experienced person.
Matlab’s help files.
The commands in Matlab have help files, which you can access from within Matlab by typing help followed by the command name, for example:
>>help plot (press the return key)
Will print the help information about the command plot. After this I will assume that you know to press the return key. The help information is the comment lines comprising the first few lines of the Matlab program. Within the program, lines containing comments in Matlab are marked by a % character. Thus you can provide help information on the programs you write in the same manner. For example, if you type help myprogram, the help lines in myprogram will be printed to the screen.
Getting out of Matlab. At the Matlab prompt, enter exit.
SOME BASIC MATLAB OPERATIONS FOR GEOPHYSICS
These exercises and example will get you started with the basic operations that are useful in geophysics. There will be some more discussion of matrices, writing to files, etc., later. You should work through this section while you are seated at a computer with Matlab running. In all cases below, you should read the help file for the functions introduced in the section.
Making a series of zeros
tz=zeros(1,100); generates a 1 by 100 matrix tz containing 100 zeros. If you don't put the ; at the end of the line, Matlab will print the array you just created to the screen. If you forget to enter one of the indexes but instead enter zeros(100), you will get an 100 by 100 array of zeros.* Try the command from the command line (the >> prompt). Verify that you have created the matrix by entering whos . Verify the values by entering:
tz. The computer will print the values.
Minihomework: also try:
tz=zeros(10,1)
whos
tz=zeros(10)
whos
Notice that the whos function tells you the names of the matrices you have created,
and also the number of dimensions and also the storage size.
Minihomework: Repeat the above example with the ones function.
Making a series of numbers representing time in seconds
Read the help for colon function. (That is, type help colon.)
For example try entering, t=0:.001:.05 Makes a one dimensional matrix, t, starting at zero, incrementing by .001, ending at .05. (51 terms).
Minihomework: Rerun the above example, but try different values for each of the three numbers (the 0, the .001 and the .05) so you understand how they function.
Minihomework: make a series of numbers representing time in seconds running from 0 to 10 seconds in increments of 0.2 seconds.
Making a series of random numbers with a "normal*" distribution.
The numbers generated by randn follow a normal or Gaussian distribution (otherwise known as a "bellshaped curve"). "Normal" means that the numbers generated follow this distribution. The reason it is called normal is that many natural processes generate this kind of distribution.
A uniform distribution may be generated by rand. See the documentation on the rand and randn for more information.
randn(1,51); creates a 1 by 51 matrix of random numbers wiuth a normal distribution. The mean of these random numbers is is 0.5.
HW: Use the command hist to plot a histogram (a probability distribution) of the values generated by rand, and comment on the shape of the plot.
% signifies that the text that follows is a comment. Contrary to popular myths and jokes, real programmers DO comment their work, and I want you to do so too. If you have to work with your program six months from now, you will wish that you wrote more comments.
INTERLUDE:
Matlab has two features that make debugging programs easier. If you type
whos, you will get a list of the current variables. This is quite handy, because your program may have created variables in a way you were not expecting. Also, typing the name of any variable
will cause the values to be listed on the screen. From this, you can conclude that Matlab keeps the variables your program has created in its workspace after you have run the program. The usefulness of this feature is that the variables are still present for you to use in another program or in the command mode. The disadvantage of this feature is that you may not want the variables there. You may want to use the same varable names and the matrix size may not be right. The clear command gets rid of all the the variables in the Matlab workspace.
Making a 50 Hz sinusoid. Enter the following lines into an .m file. You may name the file
fifty.m
t=0:.001:0.5;
f=50;
fiftyHz=sin(2*pi*f*t);
plot(t, fiftyHz)
This will generate and plot a sinusoid that is sampled every .001 seconds on the time interval 0 to 0.5 seconds. The computer will list the series of numbers you have generated if you enter
FiftyHz.
Minihomework: Modify the above short program to generate a 100 Hz cosine signal over the interval 0 to 0.05 seconds sampled every .001 second.
Minihomework 2: Remove the semicolon at the end of one of the lines of the above program. You will see why the semicolon is needed.
Notice that the sine function operates on the matrix t, term by term. Thus, Matlab saves programming effort by minimizing the need to write do loops. If you have never heard of programming with loops, you will later.
Making signal plus noise:
Add these program lines to the previous example.
r=rand(1,length(fiftyHz)); % random numbers, same number as in FiftyHz
w=1.0; % weight for the noise
mr=mean(r); % the mean of the noise
spn=fiftyHz+(w*(rmr)); % I subtracted the mean of r, to reset
% the mean to zero.
plot(spn) %plot it
In the last line, the scalar mr is subtracted from each value in the matrix r.
The weight, w, multiplies each value in the matrix rmr. The matrix fiftyHz is added to each term of the matrix (w*(rmr)). Noticice how you can mix simple arithmetic with matrix arithmetic.
Minihomework 1: Add a line to the above program to plot the function spn. Add more lines to plot r and fiftyHz. Change the value of the weight w and see the effect on the matrix spn.
Minihomework 2: Add lines to the above program to make (and plot) a second 50 Hz sinusoid plus noise with the sinusoid amplitude set to 3.0 and the noise weight set to 0.5.
Convolution.
Convolution is an operation carried out with two time series. A time series is a series of numbers (a onedimensional array or matrix), representing sampled data. Physically, for convolution, one time series can be thought of as a sampled signal and the other time series can be thought of as the impulse response of a physical system. Here is a nongeologic example of convolution: The signal might be the elevation of a car tire driven along a bumpy road, call it A. The response of the physical system might be the combined response of the car's suspension system, which is designed to absorb the bumps before delivering them to the rear end of the driver, call it B. B is the elevation of the top of the seat cushion due to a unit impulse of elevation delivered to the bottom of the car tire.
If the bottom of the car tire is subjected to the elevation time series caused by the bumpy road, then the elevation of the top of the seat cushion (the time series C) is the convolution of the signal, A, convolved with the impulse response of the suspension, B. The physical reasoning behind this is just superposition and linearity, but I will skip the details for now. The good news is that there is a Matlab function for carrying out convolution, which is called conv. Please read the documention on conv.
Please read Yilmaz pp 18 and 19 for very basic examples of convolution, autocorrelation and
crosscorrelation. For example, his computation in his Table 11 and 12 can be carried out on the Matlab command line as:
conv([1 0 0.5], [1 0.5])
ans = 1.0000 0.5000 0.5000 0.2500
Several textbooks use the following time series as an example of convolution. You could check your understanding of the conv command with these series.
4 2 1 convolved with 1 1 1 2 2 2 = 4 6 7 11 13 14 6 2
Auto and CrossCorrelation
"Seismic processing often requires measurement of the similarity or timealignment of two traces. Correlation is an operation carried out directly on the two time series that is used to make such measurements" (from Yilmaz page 18). Crosscorrelation compares one time series to another. Autocorrelation compares one time series to itself.
There is no autocorrelation routine in Matlab. To compute a autocorrelation, you must reverse* one times series and convolve it with the other time series. You reverse the time series with flipud or fliplr.
flipud means flip updown, and fliplr means flip rightleft, whichever fits the kind of matrix you have.
Here is the example of the autocorrelation of wavelet 1 from Yilmaz page 19, Table 17.
w1=[2 1 1 0 0]; % wavelet 1
conv(w1,fliplr(w1))
ans = 0 0 2 1 6 1 2 0 0
Here is an example of the crosscorrelation of two time series from Yilmaz Table 5
w1=[2 1 1 0 0]; % wavelet 1
w2=[0 0 2 1 1]; % wavelet 2
cc=conv(w1,fliplr(w2))
answer: cc= 2 1 6 1 2 0 0 0 0
You can carry out this computation either from the Matlab command line or from a little program.
Minihomework 1: Find the result if you reverse the order of crosscorrelation.
That is, change the last line to cc=conv(fliplr(w2),w1)
Answer: The output time series is reversed.
Mini homework 2: Demonstrate what happens when you reverse the order of the series used in convolution (convolution of a with b, compared to the convolution of b with a.
Answer: The output time series is not reversed.
You will discover that if the mean is nonzero, you get a triangle or trapezoid function superimposed on your output.
Minihomework: Create two simple "box car" time series , a and b as follows:
a=[ 1 1 1]
b=[1 1 1 1 1 1]
convolve them
conv(a,b),
and convolve each with itselfL:
conv(a,a)
conv(b,b).
Plot the input function and the output function. Force the horizontal axis to be the same on all the plots.
Comment on the result.
The results are memorable, and they explain the origin of the problem described above.
Minihomework: Create nontrivial signals that demonstrate the problem described above, then fix the problem by subtracting the mean from the data.
Homework: The classic simple vibroseis problem goes here. Several textbooks illustrate the principle of the use of a swept sinusoid as a seismic signal source. The reflected waveform is a composite of the swept sinusoids reflected from each subsurface layer. The transmitted wave form is correlated with the composite waveform to recover an output correlogram that resembles a conventional seismic trace. A good homwork problem is to produce all these traces in Matlab.
But, a complete statement of the problem is too lengthly to include here. An example of the traces is given in Reynolds, p 258.
Termbyterm multiply of two one by N matrices of equal length is a.*b. The operator is .* (dot star). There is also a term by term operators for divide (./) and for other operations. Note that the dot comes first.
To delay a signal, set up another array, such as 0,0,1 and convolve it with the original signal. The output will be the original delayed by two samples. That is, there will be two
zeros added to the front of the original signal. Note the syntax of the statement assigning values to a: a=[ 0 0 1] The elements of the matrix are enclosed in brackets. The values are separated by spaces, not commas.
A universal measure of signal strength is the RMS value. The initials stand for Root Mean Square. Remember that physics lab? It is computationally identical to the standard deviation. The Matlab command is std.
Odd and ends
plot(name of onedimensional matrix) plots the matrix in new plot window. The plot command works for 1 by N matrices, in which case the horizontal axis is just the index of the matrix element. There are MANY variations on the plot command. You will need to plot one set of numbers (for example fiftyHz, above on the vertical axis) versus another (like t, above, on the horizontal axis). Check the documentation with help plot.
stem(name of matrix) makes a plot with little matchsticks supporting each data point. This is sometimes used to emphasize that the points are discrete.
bar(name of matrix) produces a plot like stem, but the little matches have no heads (they are just little sticks). You can vary the width of the stick (see the help file).
subplot (m, n, l) subdivides prepares the plot window so that the next plot command plots in a subwindow. For example subplot(4,2,1) subdivides the plot window into four rows of 2 columns each, and the next plot command with plot into the 1^{st} window (the upper left hand corner). The windows are counted from left to right and from top to bottom, like frames in a Sunday comic strip.
hold on allows you to plot multiple data sets in one plot axis. Place the command before the first time you use the plot command. Then multiple plot commands will plot onto the same axis.
hold off cancels hold on and allows you to go on to another plot axis. Place this command after the last instance of plotting.
print from within Matlab prints the plot on the system printer. You can also plot by
clicking on the proper tab at the top of the plot window.
Odds and ends with almost no explanation. A few handy Matlab commands are given below. You should check the online documentation for details.
abs(absolute value of a complex number)
ang(angle of a complex number)
scale* ( [ xmin xmax ymin ymax]) scales the plot manually
title* ( adds titles and axes labels to plots)
*You enter these commands after the plot command.
PART II. MATRICES, LARGER DATA SETS, PLOTTING AND CONTOURING
Uses Matlab programs clip.m, cont4.m, plottilt.m, plotfraz.m, maglines.m, emmov.m, slicemovie.m are included on the diskette. Data sets for programs are also included.
In the first part of these notes above, I have given you very short examples, with just a little bit about working with data in matrices. In the second part of the notes, which you are now reading, I want to present methods of working with somewhat larger data sets, and do more with matrices and plotting. This part of the notes is suitable for use in a longer workshop or a termlong class.
Using Data in Matrices
The twodimensional matrix is the basic data structure of Matlab. Please review what you know about twodimensional arrays or matrices. The key to understanding matrices is the indices or the rows and the columns. My program magconts reads data into an array d. If you run the program by entering magconts, then enter whos you will see a line:
d 1700x5 68000 double array
among other lines, informing you that d is a matrix with 1700 rows and 5 columns.
You can examine this point by point if you want, entering
d(1,1) which will print the number in row 1 column 1
d(1,2) which will print the number in row 1 column 2, etc.
Usually, you will not need to refer to the individual values.
In magconts, I use data in the individual columns of d. This gives an opportunity to explain one use of the : (colon). You can use it to refer to an entire column or an entire row of a matrix.
load d.dat % load file .. numbers go into
% matrix d
x=d(:,1); % assign x to first column of d
y=d(:,2); % assign y to second column of d
z=d(:,5); % assign new z to fifth column of d
In particular, you need to know what happens when you multiply one matrix by another, and what it means to transpose a matrix. Every now and then you may need to invert a matrix to solve a set of linear equations.
The way you locate numbers in matrices is by the row and column index. The notation always refers to the row first then the column. If you have trouble remembering the order, try to remember "Roman Catholic" for Row Column. If you have trouble remembering which way the rows and columns run, try to remember that Roman buildings were held up by vertical columns.
Single dimensional matrices are a special simple kind of matrix. In Matlab, a one dimensional array or matrix can be either a matrix with one column or a matrix with one row. If you don't know which you have, it can cause problems with subsequent computations, but fortunately, it is easy to change the data from one form to the other with the transpose command.
Another very common matrix manipulation is to multiply the elements of one matrix termbyterm by the elements of another matrix or divide the elements of one matrix termby term by the elements of another matrix. This is done by preceeding the multiply or divide sign by a dot, e. g .* for term by term multipl or ./ for term by term divide.
Matlab can easily invert matrices. I will not present applications of this. If you know why it is useful to invert matrices, you can probably figure out how to do the inversion in Matlab. Two geophysical applications are:
1. finding the inverse of a wavelet for deconvolution (Yilmaz, 1987, Chapter 2).
2. estimating the minerals present from a suite of well logs e.g. Doveton (1986)
