Lab-1: Matlab/Octave and sound processing
The goal of this first lab is to gain familiarity with Matlab/Octave in the
context of sound processing. Octave is available as Free Software and can be downloaded from http://www.octave.org/download.html and its manual is available at http://www.gnu.org/software/octave/docs.html. For information on Matlab go to The MathWorks site.
[From now on the labs are presented as if you would use Octave, but there is basically no difference with Matlab]
1.1 Basic commands
The following exercises will begin your orientation in Octave.
(a) Open an "xterm" and type "octave". To be able to read and write
sound you need to download a sound from Freesound.org (ex: corto.wav). Put it in the directory where you started
(b) Explore the Octave help capability. Type each of the following lines to read about these commands:
- help plot
- help zeros
- help ones
- help -i operators
- help -i colon
(c) Use Octave as a calculator. Try the following:
- pi*pi - 10
- ans ^ 2 %<--- "ans" holds the last result
(d) Variable names can store scalars and matrices in Octave. Try the following:
- xx = sin( pi/5 );
- cos( pi/5 ) %<--- assigned to what?
- yy = sqrt( 1 - xx*xx )
(e) Complex numbers are available in Octave. Notice that the names of
some basic operations are unexpected, e.g., abs for magnitude. Try the
- zz = 3 + 4i
- exp( sqrt(-1)*pi )
- exp(j*[ pi/4 -pi/4 ])
(f) Plotting is easy in Octave, for both real and complex numbers.
The basic plot command will plot a vector yy versus a vector xx. Try
Drop the semicolons, if you want to display the values in the x,
y, and z vectors. xx.*xx denotes an element wise calculation in
contrast to a matrix multiplication xx*xx . When unsure about a
command, use help.
- x = [-3 -1 0 1 3 ];
- y = x.*x - 3*x;
- plot( x, y )
- z = x + y*sqrt(-1);
- plot( z) %<---- complex values can be plotted
1.2 Array indexing
(a) Make sure that you understand the colon notation. In particular, explain what the following Octave code will produce
- i = 2 : 4 : 17
- k = 99 : -1 : 88
- l = 2 : (1/9) : 4
- m = pi * [ 2 : (-1/9) : 0 ]
(b) Extracting and/or inserting numbers in a vector is very easy to do. Consider the following definition:
Explain the result echoed from the last three lines of the above code.
- xx = [ ones(1,4), [2:2:11], zeros(1,3) ]
(c) In the previous part, the vector xx contains 12 elements. Observe the result of the following assignment:
Now write a statement that will replace the odd-indexed elements of xx
with the constant -77 (i.e., xx(1), xx(3), etc). Use a vector indexing
and vector replacement.
1.3 Use vectors instead of loops
(a) Experiment with vectors in Octave. Think of a vector as a list of numbers. Try the following:
Explain how the last example computes the different values of
cosine without a loop. The text following the % is a comment; it may be
omitted. If you remove the semicolon at the end of the first statement,
all the elements of kset will be echoed to the screen.
- kset = -3:11;
- cos( pi*kset/4 ) %<---comment: compute cosines
(b) Vectorization is an essential programming skill in Octave. Loops
can be done in Octave, but they are not the most effcient way to get
things done. It's better to avoid loops and use the vector notation
instead. For example, the code below uses a loop to compute values of
the sine function. Rewrite this computation without using the loop (as
in the previous part).
- xx = [ ]; %<--- initialize the x vector to a null
- for k=0:7
- xx(k+1) = sin( k*pi/4 ) %-- xx(0) would fail
1.4 Script files
(a) Use an editor such as emacs or notepad (on UNIX or DOS), to
create a script file called funky.m containing the following lines:
Save the script e.g. in the folder octave_files.
- tt = -2 : 0.05 : 3;
- xx = sin( 2*pi*0.789*tt );
- plot( tt, xx ), grid on %<--- plot a sinusoid
- title('TEST PLOT of SINUSOID')
- xlabel('TIME (sec)')
(b) Run your script from Octave. To run the skript funky that you created in part (c), try
(c) Add three lines of code to your script, so that it will plot a
cosine xx2=0.5*cos( 2*pi*0.789*tt ) on top of the sine. Two graphs can
be plotted on top of each other by
funky <---will run the commands in the file
- type funky <---will type out the contents of funky.m to the screen
- which funky <---will show directory containing funky.m
"b" denotes blue color, "r" denotes red color
plot(tt, xx,"b", tt, xx2,"r").
(a) Now generate a tone (i.e., a sinusoid) and listen to
it with the sound command. The frequency of your tone should be 2 KHz
and the duration should be 1 sec. The following lines of code should be
saved in a file called mysound.m and run from the command line.
The format of the wavwrite command might be different for your version
of octave. if you get an error from the above code, try to figure out
the right syntax using "help wavwrite".
Play the sound file "sound.wav" outside Octave with any sound player,
for example "play sound.wav" from the commandline. You can use wavread
to read in a sound from a file.
- dur = 1.0;
- fs = 8000;
- tt = 0 : (1/fs) : dur;
- xx = sin( 2*pi*2000*tt );
- wavwrite( "sound.wav",xx', fs, 16)
1.5 FunctionsThe following warm-up exercises are to help you
in writing functions in Octave. Although these examples contain minor
errors, they do exemplify the correct structure and syntax for writing
(a) Find the mistake in the following function:
(b) Find the mistake in the following function:
- function xx = cosgen(f,dur)
- %COSGEN Function to generate a cosine wave
- % usage:
- % xx = cosgen(f,dur)
- % f = desired frequency
- % dur = duration of the waveform in seconds
- tt = [0:1/(20*f):dur]; % gives 20 samples per period
- yy = cos(2*pi*f*tt);
- function [sum,prod] = sumprod(x1,x2)
- %SUMPROD Function to add and multiply two complex numbers
- % usage:
- % [sum,prod] = sumprod(x1,x2)
- % x1 = a complex number
- % x2 = another complex number
- % sum = sum of x1 and x2
- % prod = product of x1 and x2
- sum = z1+z2;
- prod = z1*z2;
(c) Explain the following lines of Octave code:
- yy = ones(7,1) * rand(1,4);
- xx = randn(1,3);
- yy = xx(ones(6,1),:);
(d) Write a function that performs the same task as the following without using a for loop.
- function Z = expand(xx,ncol)
- %EXPAND Function to generate a matrix Z with identical columns
- % equal to an input vector xx
- % usage:
- % Z = expand(xx,ncol)
- % xx = the input vector containing one column for Z
- % ncol = the number of desired columns
- xx = xx(:); %-- makes the input vector x into a column vector
- Z = zeros(length(xx),ncol);
- for i=1:ncol
- Z(:,i) = xx;
(a) Explain the following lines of Octave code:
- A = randn(6,3);
- A = A .* (A>0);
(b) Write a new function that performs the same task as the
following function without using a for loop. Use the idea in part (a).
In addition, the Octave logical operators are summarized via help
- function Z = replacez(A)
- %REPLACEZ Function that replaces the negative elements
- % of a matrix with the number 77
- % usage:
- % Z = replacez(A)
- % A = input matrix whose negative elements are to
- % be replaced with 77
- [M,N] = size(A);
- for i=1:M
- for j=1:N
- if A(i,j) < 0
- Z(i,j) = 77;
- Z(i,j) = A(i,j);
1.7 Manipulating sinusoids
Generate two 3000 hertz sinusoids with different amplitudes and phases.
- x1(t) = A1 cos(2pi(3000)t + phi1)
- x2(t) = A2 cos(2pi(3000)t + phi2)
(a) Select the value of the amplitudes as follows: let A1 = 13 and
use your age for A2. For the phases, use the last two digits of your
telephone number for phi1 (in degrees), and take phi2 = -30degrees.
(b) Make a plot of both signals over a range of t that will exhibit
approximately 3 cycles. Make sure the plot starts at a negative time
so that it will include t = 0, and make sure that your have at least
20 samples per period of the wave.
(c) Verify that the phase of the two signals x1(t) and x2(t) is
correct at t = 0, and also verify that each one has the correct maximum
(d) Use subplot(3,1,1) and subplot(3,1,2) to make a three-panel
subplot that puts both of these plots on the same window. See help
(e) Create a third sinusoid as the sum: x3(t) = x1(t) + x2(t). In
Octave this amounts to summing the vectors that hold the samples of
each sinusoid. Make a plot of x3(t) over the same range of time as used
in the previous two plots. Include this as the third panel in the
window by using subplot(3,1,3).
(f) Measure the magnitude and phase of x3(t) directly from the plot. Explain how the magnitude and phase were measured.
1.8 Reading and understanding sounds
A digital sound is characterized by its sampling rate (FS) in Hertzs and the number of bits per sample (BITS).
(a) Read the sound file you downloaded in section 1.1 using the function WAVREAD.
- What is its sampling rate?
- What is the number of bits per sample?
- What is the length of the sound in seconds?
the sound waveform using PLOT, expressing the horizontal axis in
seconds and the vertical axis in normalized amplitude, from -1 to 1.
a periodic part of the sound and draw only one period. What is the
period length in seconds. What is the corresponding frequency in Hz?
and draw the amplitude envelope of the sound. The envelope should have
a value for every 300 samples of the waveform. Write an OCTAVE function.