MPTK in Python: basics

Today, I begin looking at starting to begin what I need to do in order to start and reproduce my wivigram code from yesterday using the python scientific visualization library matplotlib. First, I download the dmg, and install it. Then I read the user guide. Then I see there is something called Enthought, which provides an all-in-one working python installation, and includes matplotlib! Then I see it costs money.
Then I quickly return to the terminal.

So, in the terminal I type “python”, and away we go! To lunch.

Now that I am back, it is time for a coffee. Such is life in France.

My colleague Dan Stowell has just given me a quick tour of his new python wrapper for MPTK (available soon).
Along with this, he gives me some extremely straightforward python code to demonstrate it, which I include and comment below. First, there is the set up:

import numpy as np # load numerical python library, alias as np
from math import exp, sin, pi # load some math functions
import random # presumably some random functions
import mptk # this is Dan's new wrapper
import matplotlib.pyplot as plt # load plotting routine of matplotlib
from scikits.audiolab import Sndfile # some functions to work with outputting sound signals
from scikits.audiolab import Format
mptk.loadconfig('/usr/local/mptk/path.xml') # set up paths to MPTK

I had to install scikits.audiolab, which was really easy. Just go here, download, unpack, change to that directory, and run “sudo python install”, as described here.
Now, we are going to create a signal.

fs = 44100. # set sampling rate, and the "." means fs is floating point
basefreq = 100.
exprate = 1.2
def sweep_pos_to_freq(pos): # define a function
return basefreq * exp(exprate * (pos / fs))

The last part is a function definition, about which I learn more here.
Now we can plot something.

sigdur = 4 # this is an integer
timeposses = [pos/fs for pos in range(int(sigdur * fs))]
plt.plot(timeposses, [sweep_pos_to_freq(pos) for pos in range(int(sigdur * fs))])
plt.title("Frequency trajectory")

And this is where the Python basics comes in handy.
From this we see “range” creates a set of integers,
and so “timeposses” is a set of sampling times.
Then, we plot the function value at each of these times, add a title, and finally reveal the figure.
Now, we synthesize a signal from a sweeping sine, and plot every 10th sample.

sigmul = 0.5
noiselev = 0.05
phase = 0.
signal = [] # start with an empty vector
for pos in range(int(sigdur * fs)): # for each sample
phase += sweep_pos_to_freq(pos) / fs # integrate phase
datum = (sin(2. * pi * phase) + random.gauss(0, noiselev)) * sigmul
plt.plot(timeposses[::10], signal[::10])

We can now write our signal to a wav file:

outpath = "sweep_orig.wav"
outsf = Sndfile(outpath, "w", Format('wav'), 1, fs)

To see what that “1” means passed to Sndfile, I type “help(Sndfile)” and see it refers to the number of channels. So “outsf” is like a fileID, and we write the samples by calling the write_frames method of the instance. Apparently, we have to cast “signal” (a list) to a numpy.ndarray.
It is time in decompose!

(book, residual, decay) = mptk.decompose(signal, 'dictGabor.xml', fs, numiters=300)

This just says, “Run MP using dictGabor.xml for 300 iterations, and give me the book, residual, and energy decay.
We can now build the approximation of signal from the book.

rec = mptk.reconstruct(book)
plt.plot(timeposses[::10], rec[::10])

Now, for this chirp signal approximated by static-frequency Gabor atoms,
how do the modulation frequencies compare with the instantaneous frequency of the signal over its duration?

fig = plt.figure()
plt.plot(timeposses, np.log([sweep_pos_to_freq(pos) for pos in range(int(sigdur * fs))]), '-', color=(0.75,0.75,0.75),linewidth=3.0)
for atom in book.atoms:
timedelta = (atom['len'][0] * 0.5)/fs
xs = [atom['pos'][0]/fs - timedelta, atom['pos'][0]/fs + timedelta]
freq      = atom['freq'] * fs
freqdelta = atom['chirp'] * fs / timedelta
ys = [freq - freqdelta, freq + freqdelta]
plt.plot(xs, np.log(ys), 'k-',linewidth=2.0)
plt.ylabel('Instantaneous Frequency ($log$ Hz)')
plt.xlabel('Time (s)')
plt.title("Atoms found")

Note the space at the end of the for loop. It is essential to avoid 10 minutes of debugging and unnecessary anger.
And here is the result!
Thanks for your help Dan!


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s