In this lab, you'll measure the amplitudes and phases of the harmonics of a musical instrument. Then you'll use those harmonic amplitudes and phases to create phasors, add them together, and synthesize a musical tune.
In order to make sure everything works, you might want to go to the command line, and run
pip install -r requirements.txt
This will install the modules that are used on the autograder, including numpy, h5py, and the gradescope utilities.
First, let's load a violin note. This sample is a violin, playing a middle-C note (C5), extracted from the file Violin.arco.mf.sulG.C5G5.aiff in the University of Iowa musical instrument database (http://theremin.music.uiowa.edu/MISviolin.html).
import soundfile as sf
violin, samplerate = sf.read('violin.wav')
import numpy as np
time_axis = np.arange(len(violin))/samplerate
import matplotlib.pyplot as plt
plt.figure(figsize=(14, 5))
plt.plot(time_axis,violin)
plt.xlabel('Time (seconds)')
plt.title('Violin Playing Middle C')
Text(0.5, 1.0, 'Violin Playing Middle C')
Let's zoom in on part of it.
plt.figure(figsize=(14, 5))
plt.plot(time_axis[samplerate:samplerate+2000],violin[samplerate:samplerate+2000])
plt.xlabel('Time (seconds)')
plt.title('Violin Playing Middle C')
Text(0.5, 1.0, 'Violin Playing Middle C')
So, obviously, it's periodic. That means that it is made up of the sum of many harmonics, each with a different amplitude and phase. But how do you find its amplitudes and phases?
At this point, you should switch to Praat (https://www.fon.hum.uva.nl/praat/).
Notice that, in the above spectrum, I've selected the fourth harmonic, so you can see that it is at a frequency of 1025Hz, and has a level of 30dB. Do the same for all of the first eight harmonics, and change the following line so it correctly lists the levels of the first eight harmonics.
levels = [37.5, 14.3, 42.6, 30.3, 17, 35.8]
At this point, we'll load the file submitted.py.
The file submitted.py is the only part of your work that the autograder will see. The only purpose of this notebook is to help you debug submitted.py. Once you have revised submitted.py enough to make this notebook work, then you should go to the command line, and type python run_tests.py. Once that command returns without errors, then you can go ahead and submit your file submitted.py to the autograder. You can submit to the autograder as often as you want, but it will save you trouble if you debug as much as you can on your local machine, before you submit to the autograder.
We will use importlib in order to reload your submitted.py over and over again. That way, every time you make a modification in submitted.py, you can just re-run the corresponding block of this notebook, and it will reload submitted.py with your modified code.
Since the file is called submitted.py, python considers it to contain a module called submitted. As shown, you can read the module's docstring by printing submitted.__doc__. You can also type help(submitted) to get a lot of information about the module, including its docstring, a list of all the functions it defines, and all of their docstrings. For more about docstrings, see, for example, https://www.python.org/dev/peps/pep-0257/.
import submitted
import importlib
importlib.reload(submitted)
print(submitted.__doc__)
This is the module you'll submit to the autograder. There are several function definitions, here, that raise RuntimeErrors. You should replace each "raise RuntimeError" line with a line that performs the function specified in the function's docstring.
Now it's time for you to open submitted.py, and start editing it. You can open it in another Jupyter window by choosing "Open from Path" from the "File" menu, and then typing submitted.py.
Once you have it open, try editing the function level_to_amplitude so that its functionality matches its docstring. Here is what it's docstring says:
help(submitted.level_to_amplitude)
Help on function level_to_amplitude in module submitted:
level_to_amplitude(levels)
amplitudes = level_to_amplitudes(levels)
levels - any array or list of levels, in dB
amplitudes - a numpy array or list, of the same length, containing amplitudes
Remember that the relationship between levels and amplitudes is given by
$$\mbox{level} = 20 \log_{10}(\mbox{amplitudes})$$Your goal is to invert that equation, so that, given levels, you return amplitudes.
Once you have performed that task correctly, the following code block should print the values
[0.001 0.1 1 10 100 1000]
(it might be printed in scientific notation).
importlib.reload(submitted)
amplitudes = submitted.level_to_amplitude([-60, -20, 0, 20, 40, 60])
print(amplitudes)
[1.e-03 1.e-01 1.e+00 1.e+01 1.e+02 1.e+03]
A phasor is just a complex number, $z=ae^{j\theta}$, where $a$ is some amplitude, and $\theta$ is some phase. Open submitted.py, and modify the function create_phasors so that it does what its docstring says it should do. Here's what its docstring says it should do:
importlib.reload(submitted)
help(submitted.create_phasors)
Help on function create_phasors in module submitted:
create_phasors(amplitudes, phases)
phasors = create_phasors(amplitudes, phases)
amplitudes - an array or list of amplitudes
phases - an array or list of the same length, containing phase angles, in radians
phasors - the resulting phasors = amplitude*exp(j*phases)
If you get it working, the following block should produce the result
[ 1 1+1j 2j -2+2j ]
(it might be written in scientific notation. There might be some very small nonzero real part in the 2j term.)
importlib.reload(submitted)
phasors = submitted.create_phasors([1,np.sqrt(2),2,2*np.sqrt(2)], [0,np.pi/4,np.pi/2,3*np.pi/4])
print(phasors)
[ 1.0000000e+00+0.j 1.0000000e+00+1.j 1.2246468e-16+2.j -2.0000000e+00+2.j]
Now, let's synthesize a musical note! The input will be a list of phasors, a fundamental frequency (in Hertz), a sampling rate (in samples/second), and a duration (in seconds). The output is a musical note.
importlib.reload(submitted)
help(submitted.synthesize_note)
Help on function synthesize_note in module submitted:
synthesize_note(z, F0, Fs, d)
x = synthesize_note(z, F0, Fs, d)
z (array of length N) - an array or list of phasors, giving the amplitude and phase of each harmonic
F0 (scalar) - the fundamental frequency, in Hertz
Fs (scalar) - the sampling frequency, in samples/second
d (scalar) - the duration, in seconds
x (array of length Fs*d) - the synthesized signal
This function creates one harmonic for each element of z, then adds them together to generate x.
Remember how you use phasors:
$$x(t) = \sum_{k=0}^{N-1} \Re\left\{ z[k]e^{j2\pi (k+1) F_0 t}\right\}$$where $t$ is the time of the sample (in seconds), and $\Re\{...\}$ means "real part of ...", and the frequency of the $(k+1)^{\textrm{st}}$ harmonic is $(k+1)F_0$ (so for example, the first harmonic is at $F_0$, the second is at $2F_0$, and so on). If your code is working, the following should print one period of a sine wave, which will look something like
[ 0.00000000e+00 7.07106781e-01 1.00000000e+00 7.07106781e-01
1.22464680e-16 -7.07106781e-01 -1.00000000e+00 -7.07106781e-01]importlib.reload(submitted)
x = submitted.synthesize_note([-1j], 1, 8, 1)
print(x)
[ 0.00000000e+00 7.07106781e-01 1.00000000e+00 7.07106781e-01 1.22464680e-16 -7.07106781e-01 -1.00000000e+00 -7.07106781e-01]
The following should plot two periods of a square wave, approximated by its first nine harmonics (of which the even-numbered harmonics all have zero amplitude, and the odd-numbered harmonics have amplitude of +/-1/n and phase of zero).
importlib.reload(submitted)
x = submitted.synthesize_note([1,0,-1/3,0,1/5,0,-1/7,0,1/9],2,100,1)
plt.figure(figsize=(14, 5))
plt.plot(x)
[<matplotlib.lines.Line2D at 0x7ffe52f1e370>]
If you measured the amplitudes of your violin note accurately, up above, the following should generate an 0.5-second note that looks and sounds kind of like a violin.
Note that we're using 0 for the phases. That's because we
importlib.reload(submitted)
amplitudes = submitted.level_to_amplitude(levels)
phasors = submitted.create_phasors(amplitudes, np.zeros(len(amplitudes)))
x = submitted.synthesize_note(phasors, 255, 8000, 1)
plt.figure(figsize=(14, 5))
plt.plot(x[:500])
[<matplotlib.lines.Line2D at 0x7ffe52f87760>]
In order to make the whole note sound more like a violin, the following code multiplies the whole note by a Hanning window. That serves to give the note a gradual onset, and a gradual offset; otherwise, it would start and end suddenly.
note = x * np.hanning(len(x))
plt.figure(figsize=(14, 5))
plt.subplot(2,1,1)
plt.plot(x)
plt.title('Original note, without the window')
plt.subplot(2,1,2)
plt.plot(note)
plt.title('Windowed note, so it has some onset and offset shaping')
import soundfile
sf.write('note.wav',note/np.amax(np.abs(note)),samplerate=8000)
import IPython.display as ipd
ipd.Audio(note, rate=8000)
You've been provided with a file named note2f0.py that contains a dict, note2f0.F0. You can use it to look up the fundamental frequency of a note by name, as follows:
import note2f0
print(note2f0.F0['C4'],note2f0.F0['A4'],note2f0.F0['C5'])
261.63 440.0 523.25
Use it to create a function submitted.names_to_F0 with the following signature:
importlib.reload(submitted)
help(submitted.names_to_fundamentals)
Help on function names_to_fundamentals in module submitted:
names_to_fundamentals(names)
F0 = names_to_fundamentals(names)
names - a list of names of notes, e.g., ['D#4','G#4','F4','G4','F4']
fundamentals - a list of the corresponding fundamental frequencies, in Hertz, e.g., [311.13, 415.3, 349.23, 392.0, 349.23]
If it works, the following block should produce a result something like
[311.13, 415.3, 349.23, 392.0, 349.23]importlib.reload(submitted)
fundamentals = submitted.names_to_fundamentals(['D#4','G#4','F4','G4','F4','D#4','G#4','F4','G4'])
print(fundamentals)
[311.13, 415.3, 349.23, 392.0, 349.23, 311.13, 415.3, 349.23, 392.0]
The last part of this MP will synthesize a song. You'll do that synthesizing a sequence of violin notes, then sequencing them to create the song.
importlib.reload(submitted)
help(submitted.synthesize_song)
Help on function synthesize_song in module submitted:
synthesize_song(fundamentals, beats, bpm, Fs, phasors)
song = synthesize_song(fundamentals, beats, bpm, Fs, phasors)
fundamentals (array) - fundamental frequencies of the notes to be played
beats (array) - duration of each note, in beats, e.g., [1,3/4,1/4,1,1]
bpm (scalar) - number of beats per minute
Fs (scalar) - number of samples per second
phasors (list or array) - amplitude and phase of each harmonic
song (numpy array) - the returned song
This function should first use beats and bpm to figure out the durations of the notes, in seconds.
Then, for each note, it should:
(1) call synthesize_note, to synthesize the note
(2) call np.hanning, to create a hanning window, and multiply it by the note
(3) call np.concatenate to concatenate the note onto the song
If it works, the following blocks should synthesize the first two bars of "Hail to the Orange".
importlib.reload(submitted)
notes = ['D#4','G#4','F4','G4','F4','D#4','G#4','F4','G4']
fundamentals = submitted.names_to_fundamentals(notes)
beats = [1, 3/4, 1/4, 1, 1, 1, 3/4, 1/4, 2]
song = submitted.synthesize_song(fundamentals, beats, 96, 8000, phasors)
plt.figure(figsize=(14, 5))
plt.plot(song)
import soundfile as sf
sf.write('song.wav',song/np.amax(np.abs(song)),samplerate=8000)
import IPython.display as ipd
ipd.Audio(song, rate=8000)
If you reached this point in the notebook, then probably your code is working well, but before you run the autograder on the server, you should first run it on your own machine.
You can do that by going to a terminal, and running the following command line:
python grade.py
If everything worked perfectly, you'll get a message that looks like this:
But suppose that something didn't work well. For example, suppose you run python grade.py, and you get the following:
This error message means that the function test_synthesize_song, in the file tests/test_visible.py, failed (all of the other tests succeeded -- only that one failed).
The error message specifies that synthesize_song was called with the following parameters:
[ 220. 329.63 17.32 698.46 20.6 41.2 116.54 55. 46.25
41.2 3520. 1396.91 16.35 466.16 329.63 4186.01 82.41 246.94
311.13 311.13],
[0.1981401 0.1537068 1.28657637 1.72814063 1.34908962 1.77778565
0.20913126 0.71278756 1.08220118 1.05674697 0.81888192 1.27969944
0.76927221 1.33997218 1.38888085 0.72640863 0.70723419 1.63271861
0.95165971 1.62999141],
89,
18586,
[ 0.94385485-0.21474033j 1.05689668-0.33092966j 0.335067 +0.86765753j
0.93081609-0.15036639j 0.5054813 +0.63242576j -0.38066782+0.87194557j
0.89222872+0.6459836j -0.88226954+0.83670909j 0.10012949-0.94724238j
1.02861232-0.4692548j ]
With those inputs, it created a signal whose 4130'th sample should have been 0.16870767993167918. What value was produced instead? If you look at the line starting AssertionError, you can see that the function produced a sample with the incorrect value of 0.08435383996583959.
How can you debug this?
The first thing to do is to run your solution here, in the notebook, to see what answer you're getting.
importlib.reload(submitted)
song = submitted.synthesize_song([ 220., 329.63, 17.32, 698.46, 20.6, 41.2, 116.54, 55., 46.25,
41.2, 3520., 1396.91, 16.35, 466.16, 329.63, 4186.01, 82.41, 246.94,
311.13, 311.13],[0.1981401, 0.1537068, 1.28657637, 1.72814063, 1.34908962, 1.77778565,
0.20913126, 0.71278756, 1.08220118, 1.05674697, 0.81888192, 1.27969944,
0.76927221, 1.33997218, 1.38888085, 0.72640863, 0.70723419, 1.63271861,
0.95165971, 1.62999141],89,18586,[ 0.94385485-0.21474033j, 1.05689668-0.33092966j, 0.335067, +0.86765753j,
0.93081609-0.15036639j, 0.5054813 +0.63242576j, -0.38066782+0.87194557j,
0.89222872+0.6459836j, -0.88226954+0.83670909j, 0.10012949-0.94724238j,
1.02861232-0.4692548j ])
plt.figure(figsize=(14, 5))
plt.plot(song)
[<matplotlib.lines.Line2D at 0x7ffe4d5a9c40>]
Now let's compare that to the solution that the autograder was expecting.
The solutions the autograder was expecting are in the file solutions.hdf5. You are strongly encouraged to browse that file, in case you need to debug any problems with the autograder. Here's how you browse it:
import h5py
with h5py.File('solutions.hdf5','r') as f:
print(list(f.keys()))
['F0', 'Fs', 'amplitudes', 'beats', 'bpm', 'd', 'fundamentals', 'levels', 'note', 'phases', 'song', 'z']
As you can see, this file contains a lot of objects, created during a sample run of the solution code with random parameters. Let's plot the song that it was expecting:
with h5py.File('solutions.hdf5','r') as f:
plt.figure(figsize=(14, 5))
plt.plot(f['song'][:])
Right away, we notice that the expected solution (above) is twice the amplitude of the solution being generated by my file! Upon inspection, I discover that I accidentally multiplied all of my outputs by 0.5. Fixing that line in submitted.py fixes all of the results.
When gradescope grades your code, it will run grade.py. It will test your code using the solutions in solutions.hdf5, and using the test code tests/test_visible.py. It will also test your code using some hidden tests. The hidden tests are actually exactly the same as those in tests/test_visible.py, but with different input parameters.
You can earn up to 10% extra credit on this MP by finishing the file called extra.py, and submitting it to the autograder.
This file has just one function, called extra.fourier_series_coefficients, with the following signature:
import extra
help(extra.fourier_series_coefficients)
Help on function fourier_series_coefficients in module extra:
fourier_series_coefficients(x, F0, Fs, N)
phasors = fourier_series_coefficients(x, F0, Fs, N)
x (numpy array) - the input signal
F0 (scalar) - the fundamental frequency, in Hertz
Fs (scalar) - sampling frequency, in samples/second
N (scalar) - number of harmonics to measure
phasors (length-N array) - Fourier series coefficients
This should compute only the positive-frequency Fourier series coefficients
(for k=1 through k=N). Assume that the first sample of the input signal, x[0],
is at time 0; the second sample, x[1], is at time 1/Fs, and so on.
Instead of averaging over one period, you should average over the
whole signal length. So you should multiply x by the complex
conjugate of each harmonic, then add over the whole signal length,
then divide by the length of the signal.
Remember that the Fourier series coefficients are normally defined as
$$X_k = \frac{1}{T_0} \int_0^{T_0} x(t) e^{-j2\pi ktF_0}dt$$In order to summarize information over the entire length of the input signal, though, let's average over the entire input signal length, thus
$$X_k = \frac{1}{L} \sum_{n=0}^{L-1} x[n] e^{-j2\pi knF_0/F_s}$$where $L$ is the length of the signal, and $F_s$ is the sampling frequency (samples/second).
When you think you have it working, you can test it by running:
python grade.py
in your terminal. Yes, indeed, this is the same code you ran for the regular MP! The only difference is that, when you unzipped mp1_extra.zip, it gave you the test file tests/text_extra.py. So now, when you run grade.py, it will grade both your regular MP, and the extra credit part.
Incidentally, what is this test doing? It is testing to see if your analysis of the violin results in the right phasors. The violin is playing a middle C, which has an $F_0=261.63$. Its sampling frequency is $F_s=44100$. The results look something like this:
importlib.reload(extra)
phasors = extra.fourier_series_coefficients(violin, 261.63, 44100, 8)
print(phasors)
[ 1.12128763e-04+1.13458363e-04j -6.60467849e-07-7.31885347e-06j 6.46102923e-05+5.66118015e-05j 1.37028287e-06+4.38058235e-06j -7.61154056e-06+4.38417659e-06j -1.25282904e-06-9.66975949e-06j -2.36485716e-06+1.70836618e-06j 1.84020592e-05+1.13789076e-05j]
In order to find something a little closer to what we're used to, let's compute the levels of these eight notes. Since these are computed from the whole file, we'll add 120dB to the levels, in order to account for the effect of averaging over the entire file:
print(120+20*np.log10(np.abs(phasors)))
[44.05613714 17.32408509 38.68020224 13.23606341 18.87372823 19.78060986 9.29983876 26.70352854]
Well! These levels, averaged over the whole file, are pretty different from what we measured right at the peak using Praat. There are some similarities, though: for example, the 1st, 3rd, 6th and 8th harmonics have high amplitude, while the 2nd, 4th, 5th, and 7th harmonics are smaller.
Congratulations! That's the end of MP1. Good luck!