## David M. Koenig

Print publication date: 2014

Print ISBN-13: 9780198722908

Published to Oxford Scholarship Online: January 2015

DOI: 10.1093/acprof:oso/9780198722908.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2017. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see http://www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 17 January 2017

# (p.369) Appendix 5 A Brief Exposure to Matlab® and Octave

Source:
Spectral Analysis of Musical Sounds with Emphasis on the Piano
Publisher:
Oxford University Press

# Matlab

In the late 1980s I got my first copy of Matlab. The slim owner’s manual began with, “If you can’t bother to read this manual, start here.” In effect, Matlab was presented using “backward chaining.” The manual quickly showed you what Matlab could do and motivated you to dig into the details to figure out how you could avail yourself of some awesome computing power. As with many circa 1980 engineers, encountering Matlab was a mind-blowing experience (remember, many of us were using BASIC, Quick Basic, and Fortran on 8 MHz IBM PC 386s). The whole Matlab system in 1989 came on one $5.25 ′′$ flexible disk. Today the Matlab system consists of a massive collection of toolboxes and can cost thousands of dollars. For the sake of full disclosure, I am a Matlab author, which means that I have complimentary access to the Matlab toolboxes needed to carry out the simulations and graph constructions in this book.

Matlab m-files consist of “scripts” which you write in a BASIC-like language and which allow you to call many built-in incredibly powerful routines or functions to carry out calculations. For example, the function fft calculates the fast Fourier transform of a vector of data. Fortunately, you do not have to understand the fast Fourier transform to use the function.

The following is an example script file, MatlabExample.m, with some comments (anything after a “%” is a comment). Try reading it line-by-line to see if it makes sense. The Matlab scripts for many of the figures in the book are available on the book’s website.

• % MatlabExample.m (anything after a % is a comment)
• % This script calculates samples of a sound wave with a
• % weak fundamental frequency of 32.7 Hz and a strong
• % 2nd harmonic.
• % It plots the samples in the time domain.
• % It generates a sound via the computer speakers.
• % It calculates the line spectrum using a function that I wrote.
• % It plots the line spectrum.
• Clear % clear the memory
• close all % close all graphs
• N=2*44100; % number of sample points (all statements end with
• %              a semi-colon)
• Fs=44100; % sampling frequency in Hz
• h=1/Fs; % sampling interval in s
• t=(0:N-1)*h; % the N time samples in s (a vector):
• %              0,h,2h,3h,…,(N-1)h
• (p.370) % The above is an example of generating a vector of values with one
• % simple statement.
• f=32.7; % fundamental frequency of the wave in Hz
• f=CallInput(f,'enter fundamental frequency of the wave'); % CallInput
• %                                             is one of my scripts
• %=====generate the wave
• y=0.1*sin(2*pi*t*f)+sin(2*pi*t*2*f); % construct the wave samples
• y=y+0.001*randn(size(y)); % randn is a function that generates
• %                            normally distributed
• %                            random noise. The
• %                            argument is the size
• %                            of y
• % (y is a vector, as is t. Both have N samples.)
• % The above is another example of how a simple expression can
• % generate a vector of values.
• % ====== plot the wave in the time domain
• fh1=figure; % define the first figure handle
• set(fh1,'DefaultLineLineWidth',1.5) % specify the line width
• plot(t,y),grid % plot the wave on a grid
• title('A Wave with only Two Harmonics') % title for the graph
• xlabel('time (sec)') % x-axis label
• ylabel('amplitude') % y-axis label
• axis([0 .1 -inf inf]) % specify the extent of the axes as in
• %                  axis([xmin xmax ymin ymax])
• %                  plot t from 0 to .1 sec and put no
• %                  limits on plotting y i.e., plot from –inf
• %                  (the minimum value) to +inf (the
• %                  aximum value
• disp('press Enter to hear the sound of the wave ');pause
• %%====listen to the sound of the wave (Octave does not have soundsc)
• soundsc(y,Fs);
• disp('press Enter to calculate and plot the line spectrum
• of the wave ');pause
• % ===set the parameters for the line spectrum
• frL=0; % frequency in Hz of left-hand side of line spectrum
• frR=3*f; % frequency in Hz of right-hand side of the line spectrum
• % ===calculate the line spectrum
• % Y is the magnitude of the FFT (a vector), F is the vector of
• % frequencies
• (p.371) [Y, F]=SimpleLineSpectrum2(y,h,frL,frR); % SimpleLineSpectrum2 is one
• %                                         of my scripts
• Ymax=max(Y); % find the largest value of the line spectrum
• % plot the line spectrum
• fh2=figure; % define the second plot
• set(fh2,'DefaultLineLineWidth',1.5) % specify the line width
• plot(F,20*log10(Y/Ymax)),grid
• title('Line Spectrum of a Wave')
• xlabel('frequency (Hz)')
• ylabel('dB')
• axis([ frL frR -100 0])
• disp('press Enter to calculate and plot the cumulative line spectrum
• of the wave ');pause
• %===calculate the cumulative line spectrum
• fh3=figure;
• [Cls,fCls]=CumLineSpectrum3(y,h); % CumLineSpectrum3 is one of
• my scripts
• plot(fCls,100*Cls),grid
• title('Cumulative Line Spectrum of a Wave')
• xlabel('frequency (Hz)')
• ylabel('accumulated power in %')
• axis([frL frR 0 100])

Use the Matlab editor to save this script and give it a name, such as MatlabExample.m. Then run it by entering the name at the Matlab prompt, as in

• >> MatlabExample

To use the scripts, set up a folder, for example C:\APlaceForBookFiles, and download or copy the files on the website to it. You will need to copy SimpleLineSpectrum2.m and CumLineSpectrum3.m also. Modify your Matlab path to include this directory.

If you try to run the scripts that generate the figures you may run into trouble because many of them require large data files which are not on the website. In those cases, to get a feel for how the script works, you can substitute samples of a sine wave (or a sum of several sine waves) of your own making. If you do this, it is probably best to generate the wave in a manner similar to the following:

• N=2*44100; % number of points for 2.0 s of data sampled at 44100 Hz
• Fs=44100; % sampling frequency in Hz
• h=1/Fs; % sampling interval in s
• t=(0:N)*h; % the time samples
• f=32.7; % frequency of the wave
• y=sin(2*pi*t*f)+.5*sin(2*pi*t*2*f); % construct the wave
• %                     (a fundamental and the 2nd harmonic)

(p.372) You can save the wave as a mat-file (a Matlab data file) in C:\APlaceForBookFiles as follows:

• Save testwave y Fs

Then in the script you can read or load the wave as follows:

You will have to comment out the parts of the script that read in the original long files.

# Octave

There is a free software package called GNU Octave, or just Octave, that has some of the features of Matlab. The above example script will run under Octave, although about 3–5 times slower. The data collection software that I use to sample the pianos described in Appendix 4 runs under Octave in concert with Audacity, which is also free, but, again, much slower and with fewer features. Mat-files saved by Octave must be saved in ASCII format if they are to be read by both Matlab and Octave. At present, Octave version 3.6.4 has a problem constructing 3D graphs using the plot3 command. This means that dealing with huge files resulting from sampling a sound wave at 44,100 Hz for several seconds is difficult. For example, constructing a CLSM for 88 keys where each key has been sampled for 10 s requires significant decimation. Alternatively, one can make CLSMs for a subset of the 88 keys, say, 30 keys at a time. As further versions are developed, this problem may be solved.

Having used Matlab for more than 20 years, I am spoiled by its features and its speed. I am just starting to use Octave. I will attempt to develop a suite of Octave scripts that will do the data collection and some of the 3D graphs and that will only require the user to understand how Windows Explorer works. The reader should look on the website for the latest versions. At the time of publication the suite consists of six Octave scripts:

 CollectPianoNotes.m Collects the waves from each piano note. See the YouTube video for guidance. AnalyzePianoNotes.m Parses the waves from the wav file collected by CollectPianoNotes and stores them as ASCII files, one per note. CutOffAndPadOct.m Reads the ASCII note files, removes the low amplitude trailing part, pads to 6 s, and saves the modified note files in another directory. ConstructCLSMOct.m Generates a CLSM from the note files. Has a limit of 30 keys at a time. StudyNoteOct.m Reads a selected ASCII note file and produces a time domain plot, a line spectrum and a cumulative line spectrum. CalcSpecCenOct.m Generates a spectral centroid map for a piano.

There is a YouTube video “Collecting Piano Data Using Octave and Audacity” that describes how to use the first two scripts. Note that AnalyzePianoNotes is called by CollectPianoNotes and is therefore transparent to the user. The third script, CutOffAndPadOct, makes use of whatever (p.373) data was collected by the first two scripts. I recommend that you apply it immediately after you finish collecting your 88 notes. StudyNoteOct provides the user with an opportunity to look at the time domain plot, the line spectrum and the cumulative line spectrum of a specific note. CalcSpecCenOct provides a spectral centroid for the piano from which you collected the data. ConstructCLSMOct is a bit dicey and should be used with care. Depending on the computer, it will probably not produce a CLSM for more than 30 keys at a time. I would appreciate feedback on how well these scripts work (or not).