CIS 558 / Linguistics 525
Computer Analysis and Modeling of Biological Signals and Systems
Homework 1

Due: 2/1/2005

Things to learn:

Your answers should be in the form of a file that can be executed directly by MATLAB (or Octave). Please check that your file can actually be executed before submitting it. Email your solutions to Mark Liberman. If you want to add any commentary, use the comment character % to do so.

These exercises should work equally well in Octave and in MATLAB. If you use Octave on linux, invoke it with the "--traditional" flag to maximize MATLAB compatibility. I've tested most of this assignment with Octave 2.1.34 and MATLAB versions 4 through 6.

See this link for information about where and how to access these programs.

This set of directions do not include an adequate tutorial introduction to MATLAB, though if you are familiar with similar computer languages, you may be able to get by. For a MATLAB tutorial, you can use one of those available on-line (try google), or the information that you can get from the menu item help>>MATLAB in MATLAB 6, or the introductionary sections of the Mastering MATLAB book, which will be made available to you in Williams 623. Given the multiplicity of versions and environments, you may have some initial troubles with issues like loading data -- one of the goals of this exercise is to get you past those.

In order to execute a .m file script, load a .mat file, or read in another sort of file, the files must be in a folder or directory that MATLAB or Octave can "see". The way(s) to do this vary somewhat with the program, the release number, and the operating system. This may involve placing the file in a special directory; and/or invoking MATLAB or Octave in the directory in question; and/or using the path( ) command; and/or using the menu item file>>set path (in MATLAB 6.X). If you run into any persistent trouble with this issue, please ask for help.

In order to create a script that you can run in MATLAB (or Octave), create a text file whose name ends in ".m" (like "foo.m" or "homework1.m" or whatever), in an accessible directory or folder. Then you should be able to type

foo

in MATLAB (or Octave) and have the contents of the file executed as if you had typed them.

A. Basic 1D operations in MATLAB

The file audio1.mat contains 10,399 numbers representing about 1.3 seconds of a sound waveform sampled 8,000 times a second, assigned to a variable named S.
  1. Load this file into MATLAB or Octave. First, download the file to a place on your machine where MATLAB or Octave can find it. Do this by right-clicking on the link in your web browser, and choosing save link target as or the equivalent from the menu that appears.
    If you are running Octave on Windows, you will want to save it in the folder "octave-files" in the folder where Octave is installed -- this might be e.g. C:\Program Files\GNU Octave 2.1.36\octave_files. Then load('audio1.mat')should find it.
    If you are running Octave on Linux, then you can save it in a directory you create for the purpose (e.g. /home/myname/matlab). Then if you cd to that directory and invoke Octave as octave --traditional, the load( ) command should find it.
    If you are running MATLAB 6 on Windows, you should be able to use the import data wizard [ file>>import data... ] ; or you can add the folder where the file is stored to your MATLAB path using the path( ) command (or [ file>>set path... ] ), and then import the file using the command load('audio1').
    Other combinations of MATLAB or Octave versions and operating systems may require slightly different ways to do this.
    After you have succeeded in saving audio1.mat locally and loading it into the program, the result should be that a variable S, bound to a vector of numbers, becomes available in your workspace.
  2. Now plot the vector S. Use the function plot( ). You can execute help plot to get information on this function. You should be able just to type plot(S) and get an acceptable result.
  3. Given that S is sampled at 8 KHz, with the first sample representing time 1/8000 and the 8000th sample representing time 1, create a second vector of the samples spanning the time period from .275 to .325 seconds inclusive. Plot it. Keep in mind that putting a semicolon after a MATLAB expression will prevent its value from being printed.
  4. Create a third vector corresponding to the 10399 samples of S scaled so that the minimum value is -1 and the maximum value is 1.
  5. Calculate the serial cross-correlation of 80-sample pieces of S, starting at sample 2200, for lags between 1 and 150 samples inclusive. In other words, calculate a vector A such that

    \begin{displaymath}
\displaystyle
a_i ~=~ \frac{\sum_{j=0}^{79}{ S_{n+j} S_{n+i+...
 ... ~ \sum_{j=0}^{79} S_{n+i+j}^2}}
 ~, n~=~2200, ~1 \le i \le 150\end{displaymath}

    You can do this using the MATLAB for loop and the functions sum( ) and sqrt( ). More efficiently, try to use your knowledge of linear algebra (inner products, vector lengths)!

    As the nomenclature implies, you are calculating the correlation of a piece of S with another piece of S at a variety of offsets. Calculations of this kind are used in periodicity detection. Keep in mind that Matlab's vector indices start with 1 rather than 0.

  6. What is the maximum value of vector A? For what index value does it occur? Use help max to see how to find this out.
  7. Plot A three times, first with the X axis in samples; second with the X axis in (time lag in) seconds (i.e. samples/8000); and third with the X axis in (time lag in) Hertz (the inverse of the measure in seconds). In the last case, limit the range of values shown to 500 Hz. and below.
  8. [Extra Credit] Calculate a version of A based only on the numerator of the equation given (i.e. based only on the inner products of an 80-element stretch of S with 150 other 80-element stretches of S. Is the maximum at a different place?

B. Basic 2D operations in MATLAB

The following problems pertain to image-processing in MATLAB. An image is stored as a two-dimensional matrix, and displayed with the command image( ). The elements of the image are typically called pixels. A modified version called imagesc( ) will scale the pixel values so that the largest is mapped to white and the smallest to black.

Execute the following commands to load the image einstein.raw into MATLAB or Octave (once you have downloaded the file to a suitable local place, and ensured that its directory or folder is in your MATLAB path):


fid = fopen('einstein.raw','r')
al = fread(fid,[256,256],'uchar')';
fclose(fid)

This should leave you with a variable named al bound to a matrix containing the image data.The MATLAB 6 "import data wizard" (file>>import data) will not work with this file, since it is a binary file that does not have a format that MATLAB recognizes.

Now display the image using imagesc( ). You may also have to execute the expression colormap(gray). If the image comes out sideways, transpose the matrix using the command

al = al';
and display it again.
  1. Select a subimage of al with vertical indices in the range 70-120 and horizontal indices in the range 80-130. This operation is commonly refered to as ``cropping''. Display the subimage. Note that even though the subimage is smaller than the original, it is displayed as being the same size (Note: you can use the command size( ) to determine the actual size of the underlying pixel array). Why is this? What happens when you resize the MATLAB figure window?
  2. Create an image containing every fourth pixel (in both the horizontal and vertical directions). This operation is known as ``subsampling'' (Hint: check help colon).
  3. Create a vector containing the 25th row of al. Plot it.
  4. Compute the minimum and maximum values of al. At what locations do these occur? Plot a histogram (essentially, a probability distribution) for the pixel values of al using the Matlab function hist( ).
  5. Matlab allows us to operate on images point-by-point. Compute (and display) a new image whose pixels are the square root of the pixels of al. Look at the histogram and compare it to that of al (you may want to execute the command figure(2) to put this in another window). Compute and display an image containing a pixel of value 1 at locations where al is greater than 128, and 0 elsewhere.
  6. Load another image (call it ron) from the file reagan.raw. Compute images that are the sum and product (pointwise!) of this image and al.