From Weigel's Research and Teaching Page

A working set of notes on basic operations in MATLAB and IDL.

The target audience is students with experience with programming that are going to work on a research project.

1. Overview

1.1. IDL

• IDL for astronomers: http://www.astro.virginia.edu/class/oconnell/astr511/IDLguide.html
• IDL plotting for LaTeX [1]
• Getting started with IDL: http://www.ittvis.com/EventsTraining/RecordedWebinars.aspx
• Intro to IDL http://nstx.pppl.gov/nstx/Software/IDL/idl_intro.html
• http://www.dfanning.com/

1.2. MATLAB

• Very basic introduction: http://cds130.org/
• Getting started with MATLAB: http://www.mathworks.com/access/helpdesk/help/techdoc/
• MATLAB Vectorization: http://blinkdagger.com/matlab/matlab-tips-and-tricks-on-manipulating-1-d-arrays-vectors/
• http://www4.ncsu.edu/~mahaider/NCSU_RTG_Site/RTM_Matlab_Intro.pdf
• http://www.math.umn.edu/~lerman/math5467/matlab_adv.pdf
• http://web.mit.edu/biophysics/sbio/PDFs/R2_slides.pdf
• http://vimeo.com/15752446

2. Problems

These notes cover what is needed for basic data manipulation.

These are the types of problems that you should be able to answer quickly if you are to claim that you can program in a language:

1. Using a for loop create an array named A with values of 1, 3, 5, ..., 99.
2. Using multiple for loops, create a 10x10 matrix with all values equal to 2.
3. Using vector notation create an array with values of 1, 3, 5, ..., 99.
4. Using vector notation, create a 10x10 matrix with all values equal to 2.
5. Using a for loop and an if statement, modify the array from problem 1. so that all values above 50 are set to zero. Call this array Am.
6. Using a for loop, count the number of elements in Am that have the value of zero.
7. Using a for loop, compute the mean and standard deviation of Am.
8. Repeat problems 1.-7. except don't use any for or while loops. (That is, use vector notation and functions such as mean, std, find.)
9. Use the function randn to create an array named R of length 100.
10. Plot a histogram of R with 10 bins using the function hist.
11. Plot a histogram of R with 20 bins using the function hist.
12. Plot a histogram of R with bins of width 1 centered on -3.5, -2.5, ..., 2.5, 3.5.
13. Create an array named Rm which only has values from R that are above 0.0.
14. Use the two outputs of the function sort to create a matrix with two columns. The first column should be values from R arranged in ascending order. The second column should be the location of the value in the array R. (For example, if R had values of 0.1, -0.1, then the matrix would have row 1 with the values -0.1, 2 and row 2 with values of 0.1, 1.)

3. Executing Scripts

4. Syntax Notes

• MATLAB code is case sensitive. A variable named Variable is not the same as variable and a function named MyFunction is not the same as MyFunction.
• Although, you may give a variable the same name as a function in MATLAB, it is not a good idea.
• IDL code is not case sensitive.

5. Plotting

6. Movies

7. Saving Images

7.1. MATLAB

7.1.1. Basic

To make your plots readable in a presentation, use

xlabel('m','FontSize',16);
ylabel('b','FontSize',16);
set(gca,'FontSize',16);
print -r600 -dpng filename.png

7.1.2. Advanced

Saving png or jpeg images in MATLAB often leads to fuzzy images because of the default settings they use in their conversion routines. In MATLAB, use

print -depsc filename.eps

and then convert the eps file with the ImageMagick convert utility that is installed by default in many Linux operating systems if the png image is not to your liking. By saving an eps file, you always have the vector eps file and can later re-rasterize it to a different size. I use the following function as a replacement to print -deps filename.eps:

function plotcmds(tmp)

 set(0,'defaultaxesfontsize',16); 
 set(0,'defaulttextfontsize',16);
 
 fprintf('plotcmds: plotting %s.png and .eps\n',tmp);
 eval(sprintf('print -depsc %s.eps',tmp))
 system(sprintf('convert -quality 100 -density 100 %s.eps %s.png',tmp,tmp));

If you do not have convert you may just want to use "print screen" to rasterize an eps file. Also, it is better to use gs than convert because convert is a wrapper to gs for many input/output format conversion pairs, and like MATLAB, convert sometimes uses poorly-chosen and difficult-to-modify settings when calling gs. For more information, see [2].

Note that what you see on the screen will often differ from what you see in the postscript file. From help print:

A note on renderers: when printing figures, MATLAB does not always use the same renderer as on screen. This is for efficiency reasons. There are cases, however, where the printed output is not exactly like the screen representation because of this. In these instances specifying -zbuffer or -opengl will more likely give you output that emulates the screen.

7.2. IDL

See also http://www.dfanning.com/graphics_tips/weboutput.html

filename = 'filename'
!P.FONT = 0
!P.REGION = [0.10,0.10,0.90,0.90]
set_plot,'PS'
device,filename=filename+'.eps', encapsulated=1, /helvetica,font_size=16,/bold
 
xdata = [0,1,2]
ydata = [0,1,2]
 
plot,xdata,ydata,psym=7,                       $
 xrange=[0,10],yrange=[0,20],/xstyle,/ystyle,  $
 xtitle='xlabel',ytitle='ylabel',title='title' 
 
device, /close
set_plot,'X'
 
command='convert -quality 100 -density 100 '+filename+'.eps '+filename+'.png'
spawn,command
print,'wrote '+filename+'.eps and '+filename+'.png'
 
stop
end

8. Basic Input/Output

8.1. Output

Writing a space-delimited text file:

>>data = [1 3;2 4];
>>save -ascii data.matlab.txt data
>>!head data.matlab.txt

IDL>data = [[1,2],[3,4]]
IDL>openW, 101, 'data.idl.dat'
IDL>printf, 101, FORMAT='(%"%0.7e %0.7e")',transpose(data)
IDL>close, /all
IDL>$head data.idl.dat

8.2. Input

See also #DataFiles for complex examples.

Reading a space or comma-delimited text file:

>>load data.txt;
>>data
data =

    1     3
    2     4
 
>>whos data
 Name      Size            Bytes  Class     Attributes

 data      2x2                32  double              

Notice that MATLAB prints the contents of the matrix data to the screen as integers, but whos reports that the values are doubles.

IDL>OPENR, lun, 'data.dat', /GET_LUN
IDL>data = FLTARR(2, 2)
IDL>READF, lun, data
IDL>close, lun
IDL>print,data

In IDL you tell READF about the structure of the file (two rows by two columns) by declaring a variable with the appropriate size. There are easier ways to read ASCII files in IDL. See

• http://www.dfanning.com/tips/ascii_column_data.html
• http://www.msi.umn.edu/software/idl/tutorial/idl-rwd.html#Reading

9. Vectors

The syntax for creating a vector is

>> A = [1,2,3,4]

IDL> A = [1,2,3,4]

Note that the commas are optional in MATLAB, but it is a good habit to include them as they are required in many other programming languages. To view A again, just type "A" on the command line and hit enter.

>> A

Put an ";" at the end of the line and hit enter. The semicolon suppresses the echoing of the input. In long scripts, you will want to terminate commands with a semicolon except in places where you want to do debugging.

IDL>print,A

Press the up arrow key to repeat the command where you defined A. To view element 2, type

>> A(2)

IDL>print,A(2) ; 

or

IDL>print,A[2] ; preferred way according to documentation

In many languages (C, Fortran, IDL) you would enter a ".0" following the integers when you created A if you meant for the elements to be floating point. That is, you would have entered

>> A = [1.0,2.0,3.0,4.0];

so that A(1)/A(2) = 0.5 and not zero (to verify this, try entering "print, 1/2" in IDL or printf("%f",1/2) in a C program. In Matlab, the result of 1/2 is 0.5, so the ".0" is not needed.

Question: What is a shorthand way to create the vector

B = [1,2,3,4,1,2,3,4,1,2,3,4]

Answer:

 >> B = [A,A,A];
IDL> B = [A,A,A]

10. Matrices (2-D)

To create a matrix, type

>> A = [11,22;33,44];

Note that here the semicolon has a different meaning than "suppress echoing of output": it is used to indicate a new row in a matrix.

IDL>A = [[11,22],[33,44]]

To view A, type

>> A

A =

    11    22
    33    44

or

IDL>print,A
     11      22
     33      44

To access the "44" using subscript notation

 >> A(2,2)

ans =

    44

 IDL> print,A(1,1)
       44

and using array notation,

 >> A(4)

 ans =

    44

 IDL> print,A[3]

MATLAB uses column major ordering, which means

>> A(2)

returns 33 and not 22. For convenience, you may use the functions SUB2IND compute the integer (i.e., 4) that corresponds to a given subscript (i.e., 2,2). IDL uses row major ordering so the zeroth element is 11, the second element is 22, etc.

IDL>print,A[1]

returns 22 and not 33.

Question: Use a vector to access the 11 and 44 (the diagonal elements) from A.

Answer

 >> A([1,4])
 IDL> A[[0,1]]]

or

 >> [A(1,1),A(2,2)] 
 IDL> [A[0,0],A[1,1]]

Question: How would you create this matrix?

11 22 11 22
33 44 33 44

Answer:

>> [A, A]

Question: How would you create this matrix?

11 22
33 44
11 22
33 44

Answer:

>> [A ; A];

11. Functions

To find out about a function, type

>> help function_name

For example

>> help size

To call a function, simply type

>> size(A)

or

>> B = size(A)

In the previous section I said that you could may access the diagonal elements of A using two forms. To verify that A(2,2) is the same as A(4), use SUB2IND, which requires information about the size of A and the requested element.

>> help sub2ind
>> sub2ind(size(A),2,2)

Question: Use the function FLIPLR to create the vector

B = [1,2,3,4,4,3,2,1];

using

A = [1,2,3,4];

12. Operators

12.1. The magic colon

• Context 1: To build a vector

>> A = [1:4];

gives

A = [1,2,3,4];

• Context 2: selecting all elements of a row or column. If

>> B = [A ; A];

select the first column of B by entering

>> B(:,1)

To select the first row of B, type

>> B(1,:)

• Context 3: to create a regular array

>> A = [1:2:5]

gives

A = [1,3,5]

So that [a:step:b] means create the array [a, a+step, a+2*step, ..., b];

Question: What is the last number created by this command?

>> A = [1:2:4]

Answer

3

So the last number created by the magic colon in this context is always less than or equal to b.

Question: Use the magic colon to create the vector

[1,-1,-3,-5]

Answer

[1:-2:-5]

Question: Use the magic colon to access every other row of the matrix A created by

>> A = rand(5,5);

Answer:

A([1:2:5],:)

• Context 4: To flatten a matrix

Suppose

>> A = rand(5,5)

To turn A into a vector, type

>> A(:)

Question: What are the dimensions of A(:)?

Answser:

size(A) = [25, 1]

12.2. Addition and Subtraction

• Scalar addition

If

>> A = [1, 3 ; 2, 4];

then

A(1) + A(2) = 3

• Matrix addition

A + A = (1+1)  (3+3) = 2 6
        (2+2)  (4+4)   4 8

Matrix addition of A and B requires that size(A) == size(B).

12.3. Multiplication and Division

• Scalar multiplication

If

A = [1, 2 ; 3, 4];

then

A(1)*A(3) =  3

or using the alternative addressing method,

A(1,1)*A(1,2) = 3

• To do element-by-element multiplication, use ".*". The "." indicates that the multiplication should be element-by-element and not regular matrix multiplication.

Question: If

A = 1 3
    2 4

what would you expect

A.*A to be?

Answer:

1 9
4 16

Question: What is A*A? Answer: It is the matrix product.

 7 15 
10 22

Question: What would you expect

A./A

to be? Answer:

1 1
1 1

Note this matrix can be created using the function ONES.

13. Looping

for i=[1:4]
   A(i) = i;
 end

Creates

A = [1,2,3,4];

Note that the square braces in "i=[1:4]" are optional and are usually omitted in this context. The square brace in statements such as

>> A = [1:4];

is also optional, but is usually not omitted in this context.

Question: What is A that results from the following loop?

for i = 3:3:8
  A(i) = i;
end

Answer:

A = [3,6];

Note that as discussed for the magic colon, [3:3:8] gets expanded to [3,6] not [3,6,9].

14. Vectorization

You have probably heard about vectorization. In scripting languages, vectorization generally leads to decreased execution time. Vectorization _usually_ leads to code that is easier to read. For example,

for i = [1:4]
  A(i) = i;
end

gives the same result as

A = [1:4];

or

A([1:4]) = [1:4];

Question: Vectorize the following:

for i=1:100
   A(i) = 101-i;
end

Answer:

A([1:100]) = [100:-1:1];

or

A = [100:-1:1];

Question: Given

>> A = reshape([1:100],10,10);

write a for loop that creates a 1-D array that contains the diagonal elements of A. Do the same without a for loop.

Question: Create A in the following two code blocks without using a for loop. Hint: You may need to use one or more of the following functions: ONES, RESHAPE, REPMAT, CUMSUM. You may find this tech note useful: html

for j=1:5
for i=1:5
  A(i,j) = i;
end
end

A solution:

Av = repmat([1:5]',1,5)
sum(Av(:)-A(:)) % Should be zero

for j=1:5
for i=1:5
  A(i,j) = i+j;
end
end

A solution

Av = cumsum([ [2:6] ; ones(4,5) ])
sum(Av(:)-A(:)) % Should be zero
% or
Av = repmat([1:5]',1,5) + repmat([1:5],5,1)
sum(Av(:)-A(:)) % Should be zero
% or (Oscar found this)
Av=hankel(2:5,5:10)
sum(Av(:)-A(:)) % Should be zero

A collection of tricks for reading data files and pre-processing data, based on issues that students encountered.

15. Data Files

15.1. ASCII Files With Non-numeric characters

MATLAB can easily read ASCII files with values that are separated by spaces and tabs and records that are separated by newlines. For example, a file named a.txt with contents

0 1 2
3 4 5

can be read using (the data will be placed in a matrix named a):

load a
a

Often data files have additional characters, for example, : and /. Such files cannot be read directly using load.

8/11/2012	14:00:00.014.262	158.866	80.3515	223.167
8/11/2012	14:00:00.029.887	158.885	80.285	223.215
8/11/2012	14:00:00.045.512	158.948	80.1195	222.661
8/11/2012	14:00:00.061.137	159.011	80.3179	222.358
8/11/2012	14:00:00.076.762	158.837	80.3302	222.826
8/11/2012	14:00:00.092.387	159.014	80.3052	223.354
8/11/2012	14:00:00.108.012	158.999	80.226	223.069
8/11/2012	14:00:00.123.637	159.055	80.1786	222.017
8/11/2012	14:00:00.139.262	158.909	80.3855	222.247
8/11/2012	14:00:00.154.887	158.833	80.3708	223.342

In the following, the entire file with one call. The file could be read line-by-line, but this will be slower.

fid = fopen('a.txt');
s = fscanf(fid,'%c'); % Read entire file at once.
fclose(fid);
s = regexprep(s,'/',' '); % Replace all slashes with a space
s = regexprep(s,':',' '); % Replace all colons with a space
% Replace patterns such as 11.98.987 with 11 98 987
% Use a program such as http://www.regexe.com/ to figure out regular expression and replacement pattern
s = regexprep(s,' ([0-9][0-9])\.([0-9][0-9][0-9])\.([0-9][0-9][0-9])',' $1 $2 $3');
% string now contains only spaces and newlines.  Convert to array using str2num.
d = str2num(s);

15.2. Reading Many Files

If the data you are interested in span many files, the following examples contain some techniques that may be used. First, create a script named getdata.m that downloads and saves the files of interest.

The following example reads the contents of the files in subdirectories of ftp://aftp.cmdl.noaa.gov/data/ozwv/SurfaceOzone/BAO/ that start with bao_o3_6m_hourlymean.

getdata.m

urlo = 'ftp://aftp.cmdl.noaa.gov/data/ozwv/SurfaceOzone/BAO/';

for yr = 2008:2009
    if yr == 2008
       m = 7; % In 2008, first file is July
    else
       m = 1;
    end
    for m = 1:12
        fname = sprintf('bao_o3_6m_hourlymean_%02d_%d.dat',m,yr); % Form the filename
        url = sprintf('%s/%d/%s',urlo,yr,fname); % Append filename to base url
        if ~exist(fname) % If file not already downloaded
            fprintf('Attempting to download %s\n',url);
            [f,stat] = urlwrite(url,fname); % Save file in current directory  
            if (stat == 1)
                fprintf('Wrote %s\n',fname);
            end
            if (stat == 0)
                fprintf('Failure\n');
            end   
        end
        fclose(fid);
    end
end

The files have the following structure (full file). In this case, each line is read.

                                                    MIXING RATIO OF SURFACE OZONE IN PARTS PER BILLION
                                                    STATION BAO  6m AGL
                                                    MONTH    JUL  YEAR 2008
                                                    EQUIPMENT  TEI O3 ANALYZER MODEL 49I

   GMT     1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19    20    21    22    23    24  MEAN   MAX
  DATE
   1   999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9
   2   999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9
   3   999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9
   4   999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9 999.9

readdata.m

clear;
i = 1;
linea = '';
for yr = 2008:2009 % Loop through years
    if yr == 2008
       mo = 7;
    else
       mo = 1;
    end
    
    for m = mo:12 % Loop through month
        fname = sprintf('bao_o3_6m_hourlymean_%02d_%d.dat',m,yr);

        fid = fopen(fname); % Open file in current directory
        for k = 1:7 % Read first seven lines and ignore.
            line = fgetl(fid);
        end
        while 1
            line = fgetl(fid); % Read a line.
            if ~ischar(line), break, end % Reached end of file.
            if (isempty(regexp(line,'[A-Z]'))) % Line does not have A-Z or a-z
                line(6:end-length(' 999.9 999.9'))
                break
                % Keep all parts of line after day of month
                linea = [line(6:end-length(' 999.9 999.9')),linea];   % Store the line and convert from string to decimal.
                i = i+1;
            end
        end
        fclose(fid);
    end
    
end

d = str2num(linea);
h = [1:length(d)];
d(d==999.9) = NaN;
plot(h,d);
xlabel('Hours since 08/01/2008');

15.3. Filling Gaps

The following data set has a gap at t=5.

t = [0,1,2,3,4,6,7];
x = [0,1,0,1,3,4,5];

For a long array, the maximum and minimum gap size can be found using

max(diff(t))
min(diff(t))

The location of the first gap can be found using

I = find(diff(t)>1);
I(1)

and to use method (a), compute the raw periodogram using

tr = t(1:I(1))
xr = x(1:I(1))

To use method (b), perform linear interpolation using the interp1 function to compute interpolated values from time t(1) to t(end):

yi = interp1(t,x,[t(1):t(end)])

16. IDL/MATLAB translation

Note - many commands not tested (especially IDL).

file = 'data.dat'
nc   = 3;
nr   = file_lines()
row  = findgen(nc,nr)
openr, lun, file, /get_lun
data = fix(fltarr(nc,row[fix(nr, type=2)]))
readf, lun, data
close, lun

% functionname.m
function [a,b] = functionname(x,y)

if (nargin < 2)
  y = 1;  
end

a = x+y;

if (nargout > 1)
 b = subfunctionname(x)
end

function c = subfunctionname(x)
c = x+1;

;functionname.pro
function functionname, x, b=b

  if (NUM_PARAMS < 2)
    y = 1;  
  end

 if keyword_set(keystring) then b = 1
 return, x+1
end 

for k =1:10,x(k) = k;,end
% or
for k =1:10
  x(k) = k;
end

for k=1,10 do x[k]=k
; or
for k=1,10 do begin
  x[k] = k
endfor

if (1 == 1)
  A(1) = 1;
elseif (1 == 2)
  A(1) = 2;
else
  A(1) = 3;	
end

if (1 EQ 1) then begin
  A[0] = 1
endif else begin
  A[0] = 2
endelse

 switch lower(method)
   case {'method1','method2'}
     disp(method);
   case 'method3'
     disp(method);
   case 'method4'
     disp(method);
   otherwise
     disp(method);
 end

case method of
  'method1': print, method
  'method2': print, method
  'method3': print, method
  else: print, method
endcase

17. Performance Comparisons

17.1. IDL/MATLAB FFT Compare

Conclusion: in IDL, select your array length carefully (code). See [5].