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 setup.py 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") plt.show()
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 signal.append(datum) plt.figure() plt.plot(timeposses[::10], signal[::10]) plt.title("Signal") plt.show()
We can now write our signal to a wav file:
outpath = "sweep_orig.wav" outsf = Sndfile(outpath, "w", Format('wav'), 1, fs) outsf.write_frames(np.array(signal)) outsf.close()
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.figure() plt.plot(timeposses[::10], rec[::10]) plt.title("Reconstructed") plt.show()
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.hold(True) 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.5)/fs xs = [atom['pos']/fs - timedelta, atom['pos']/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.grid(True) plt.xlim((0,sigdur)) plt.ylabel('Instantaneous Frequency ($log$ Hz)') plt.xlabel('Time (s)') plt.title("Atoms found") plt.hold(False) fig.savefig('test.eps') plt.show()
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!