Introduction version 04/06/03
|
|  Yüklə 141.5 Kb.
|
INTRODUCTION version 04/06/03Matlab 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 MATLABThis 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.
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 tz. The computer will print the values. Mini-homework: 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. Mini-homework: Repeat the above example with the ones function.
Making a series of numbers representing time in secondsRead 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).
Mini-homework: 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 "bell-shaped 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. Mini-homework: Modify the above short program to generate a 100 Hz cosine signal over the interval 0 to 0.05 seconds sampled every .001 second. Mini-homework 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*(r-mr)); % 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 r-mr. The matrix fiftyHz is added to each term of the matrix (w*(r-mr)). Noticice how you can mix simple arithmetic with matrix arithmetic. Mini-homework 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. Mini-homework 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 one-dimensional 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 non-geologic 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 1-1 and 1-2 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 Cross-Correlation "Seismic processing often requires measurement of the similarity or time-alignment 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). Cross-correlation 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 up-down, and fliplr means flip right-left, whichever fits the kind of matrix you have. Here is the example of the autocorrelation of wavelet 1 from Yilmaz page 19, Table 1-7. 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. Mini-homework 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.
Answer: The output time series is not reversed. You will discover that if the mean is non-zero, you get a triangle or trapezoid function superimposed on your output. Mini-homework: 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. Mini-homework: Create non-trivial 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. Term-by-term 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 endsplot(name of one-dimensional 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 1st 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 on-line 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 CONTOURINGUses 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 term-long class. Using Data in Matrices The two-dimensional matrix is the basic data structure of Matlab. Please review what you know about two-dimensional 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.
% 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 term-by-term by the elements of another matrix or divide the elements of one matrix term-by- 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) |