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 octave.

(b) Explore the Octave help capability. Type each of the following lines to read about these commands:

help
help plot
help zeros
help ones
help -i operators
help -i colon

pi*pi - 10
sin(pi/4)
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 )
ans

(e) Complex numbers are available in Octave. Notice that the names of some basic operations are unexpected, e.g., abs for magnitude. Try the following:

zz = 3 + 4i
conj(zz)
abs(zz)
angle(zz)
real(zz)
imag(zz)
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 the following:

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

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.

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:

xx = [ ones(1,4), [2:2:11], zeros(1,3) ]
xx(3:7)
length(xx)
xx(2:2:length(xx))

Explain the result echoed from the last three lines of the above code.

xx(3:7) = pi*(1:5)

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:

kset = -3:11;
kset
cos( pi*kset/4 ) %<---comment: compute cosines

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.

(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
end
xx

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:

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)')

Save the script e.g. in the folder octave_files.

(b) Run your script from Octave. To run the skript funky that you created in part (c), try

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

(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

plot(tt, xx,"b", tt, xx2,"r").

"b" denotes blue color, "r" denotes red color

1.4 Waves

(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.

dur = 1.0;
fs = 8000;
tt = 0 : (1/fs) : dur;
xx = sin( 2*pi*2000*tt );
wavwrite( "sound.wav",xx', fs, 16)

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.

1.5 Functions

The 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 functions.

(a) 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);
end

(b) Find the mistake in the following function:

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;
end

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;
end

1.6 Vectorization

(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 relop.

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;
else
Z(i,j) = A(i,j);
end
end
end

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 amplitude.

(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 subplot.

(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?
Draw the sound waveform using PLOT, expressing the horizontal axis in seconds and the vertical axis in normalized amplitude, from -1 to 1.
Find 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?
Compute 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.