MP1: Phasors

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.

---

Part 1: Measuring harmonic spectral levels (in dB)

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

We can also listen to it using IPython.display:

import IPython.display
IPython.display.Audio(violin,rate=samplerate)

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/).

1. In the "Praat Objects" window, click the "Open" button, and choose "Read from file." In the file browser that pops up, navigate to the directory containing the file "violin.wav," and choose it. The words "Sound violin" should show up in your "Praat Objects" window.
2. Click on "Sound violin," then click "View & Edit".
3. Click your mouse near the peak of the waveform, then click on the "Spectrum" button, and choose "View Spectral Slice". You should see something that looks like the following.

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]

---

Part 2: Converting levels to amplitudes

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]

---

Part 3: Creating phasors

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+0j  1+1j  0+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]

---

Part 4: Synthesizing a musical note

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 0x7f8a367ea0b0>]

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, so the waveform will look like it has the same approximate structure as a violin waveform, but with different symmetry -- the real violin is asymmetric in the time domain, but with 0-phase, the result has even symmetry in the time domain. Even symmetry means that $x(t)=-x(-t)$, just like a cosine.

The human ear is not very sensitive to phase. If you listen to the waveform, it should sound a lot like a violin.

importlib.reload(submitted)
print('The levels are',levels)
amplitudes = submitted.level_to_amplitude(levels)
print('The corresponding signal amplitudes are',amplitudes)
phasors = submitted.create_phasors(amplitudes, np.zeros(len(amplitudes)))
print('The phasors are',phasors)
x = submitted.synthesize_note(phasors, 255, 8000, 1)
plt.figure(figsize=(14, 5))
plt.plot(x[:500])

The levels are [37.5, 14.3, 42.6, 30.3, 17, 35.8]
The corresponding signal amplitudes are [ 74.98942093   5.18800039 134.89628826  32.73406949   7.07945784
  61.65950019]
The phasors are [ 74.98942093+0.j   5.18800039+0.j 134.89628826+0.j  32.73406949+0.j
   7.07945784+0.j  61.65950019+0.j]

[<matplotlib.lines.Line2D at 0x7f8a36d9ff40>]

import IPython.display
IPython.display.Audio(x,rate=8000)

In order to make the whole note sound more like a violin, the following code multiplies the whole note by a Hamming window. That serves to give the note a gradual onset, and a gradual offset; otherwise, it would start and end suddenly.

note = x * np.hamming(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)

---

Part 5: Converting note names into fundamental frequencies

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]

---

Part 6: Synthesizing a Tune

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)

---

Part 8: Grade your code on your own machine before submitting it!

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:

..... ---------------------------------------------------------------------- Ran 5 tests in 0.275s OK

But suppose that something didn't work well. For example, suppose you run python grade.py, and you get the following:

....F ====================================================================== FAIL: test_synthesize_song (test_visible.TestStep) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/jhasegaw/Dropbox/mark/teaching/ece401/ece401labs/21fall/mp1/tests/test_visible.py", line 71, in test_synthesize_song self.assertAlmostEqual( AssertionError: 0.16870767993167918 != 0.08435383996583959 within 7 places (0.08435383996583959 difference) : *** 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 ]), sample 4130, should be 0.16870767993167918 ---------------------------------------------------------------------- Ran 5 tests in 0.226s FAILED (failures=1)

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 0x7f8a489ec1c0>]

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. First, make sure that you have the module h5py installed on your computer:

!pip install h5py

Collecting h5py
  Downloading h5py-3.9.0-cp310-cp310-macosx_10_9_x86_64.whl (3.2 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.2/3.2 MB 12.0 MB/s eta 0:00:0000:0100:01
Requirement already satisfied: numpy>=1.17.3 in /Users/jhasegaw/anaconda3/lib/python3.10/site-packages (from h5py) (1.24.3)
Installing collected packages: h5py
Successfully installed h5py-3.9.0

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'][:])

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.

---

Extra Credit

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!