Halloween Raspberry Pi project

Halloween is a great time to do projects if you have kids.  There is a lot of fun to be had building animated things for displays or costumes.  In our case, Madeline and I collaborated on a great costume for our son, who wanted to go as a Robot.  She found some shiny material and sewed it up, and I build an animated control panel.

The Rπ0 controls a IIC-connected OLED display (128x32) and 9 LEDs via the GPIOs. I could theoretically have used an Arduino - but the OLED display is 3.3V only, and I only have 5V Arduinos (I know I can just underclock them, but not during development).  Thickness was also a factor - the Rπ0 is VERY flat.

The panel is built mostly of Moosgummi (craft foam rubber), and I mounted my Raspberry Pi 0 into it.  The Rπ0 and all the other stuff is mounted in a sandwich of Moosgummi sheets with strips of Moosgummi around the edge.  Power and development (eg. changing the scrolling message or the blink pattern) is done trough a gap, through which a (short) USB cable is fed.  The LEDs are connected trough 150 Ohm resistors directly to the GPIOs.

The SD card has Raspbian Stretch Lite on it, to which I added the Adafruit SSD1306 Python library.  The Python script (placed in a Gist here) running the display and the LEDs sits in the (FAT-formatted) boot partition, and is called from rc.local.  Once that was all in place, I used the Read-Only Raspbian script to make the system read-only; this way I can just pull the power once the trick-or-treating is done.  The script on the boot partition can still be edited at will by remounting /boot read-write (or pulling the SD card and editing it on a different computer!)  The Rπ0 acts as a USB Ethernet gadget, and I can just ssh into it.

Life lesson learned: Mounting blinking LEDs on your kids makes them very easy to find at night! Note to self: next time, add lights to the back (of the kids) as well.

NNResample package on PyPI

In the previous post, I described my filter generation method for scipy.signal.resample_poly, where the filter is designed such that the first null if the filter is on Nyquist.  For easier inclusion into my own projects, I made it into a package as a self-standing function, and then decided to make that my first submission to PyPI.  The package depends only on numpy and scipy>=0.18.0.

So, if you want to use my resampler, you can now simply do `pip install nnresample`. To report bugs, the link to the repository is https://github.com/jthiem/nnresample.

I've included one "special feature": since it takes a little time to generate the filter (to search for the null and regenerate the filter afterwards) the generated filter is cached in a global for re-use.

As for the name, since I'm calling it "Null-on-Nyquist Resampler", the obvious name should be non-resample, but I think that could be confusing...

A library of Gammatone Filterbanks in Python

I have already programmed gammatone filterbanks several times in MATLAB, but in for some recent work I needed a specific one (the One-Zero Gammatone filterbank) [Katsiamis, 2006] in Python.

Gammatone impulse responses, in time domain

In addition, I wanted my "old" FIR GTFB (from my Ph.D. Thesis) to be in there as well as the GTFB of a former colleague here in Oldenburg [Chen, 2015].  So I packed them all up in a library - I also wanted to learn how to make a PyPI compatible Python library.

This is a work in progress. The gammatone filterbanks mentioned are there, but momentarily I don't have time to polish the code; furthermore, I'd like to eventually add the Slaney GTFB and the Hohmann GTFB - not to mention the structure needs more cleanup.  Consider this v0.01, for educational purposes. You can find it on github.

Audio Resampling in Python

 (Note: this post was written as a Jupyter Notebook which can be found with the Python code at https://nbviewer.jupyter.org/gist/jthiem/0bb8b2296c8dbd0b83bfab8270a8eb24.)

UPDATE: You can resample with nothing but scipy and one extra function now. See this new post.
I now have written a package that you can install using pip.

One of the most basic tasks in audio processing anyone would need to do is resampling audio files; seems like the data you want to process is never sampled in the rate you want. 44.1k? 16k? 8k? Those are the common ones; there are some really odd ones out there.

Resampling is actually quite hard to get right. You need to properly choose your antialias filters, and write a interpolation/decimation procedure that won't introduce too much noise. There have been books written on this topic. For the most part, there is no single univeral way to to this right, since context matters.

For a more detailed discussion on sample rate conversion see the aptly named "Digital Audio Resampling Homepage" [https://ccrma.stanford.edu/~jos/resample/]. Then look at [http://src.infinitewave.ca/] for a super informative comparison of a ton of resampling implementations. No seriously, go there. (It inspired me to compare the methods below using a sine sweep.)

MATLAB has the built-in resample function. Everyone uses it. It works, but also sucks in many aspects - a much better function is in the DSP toolbox, but as near as I can tell, noone uses it (even when copious other parts of the DSP toolbox are being used!)

But I'm not talking about MATLAB, I'm talking about Python here, and it doesn't have a default resampling method - except if you're doing any kind of numerical computing, you'll be using numpy and probably scipy --- and the latter has a built-in resample function.

To get ahead of myself a bit, scipy.signal.resample sucks for audio resampling. That becomes apparent quite quickly - it works in frequency domain, by basically truncation or zero-padding the signal in the frequency domain. This is quite ugly in time domain (especially since it assumes the signal to be circular).

There are two alternatives I want to point at that are "turn-key", that is, you just have to specify your resampling ratio, and the library does the rest. I know of a few others, but at the time of writing this they didn't seem as easy to use (eg. just doing the interpolation/decimation, but not calculating a filter.) If there are other notable libraries I should know of, please mention it in the comments of the blog version of this notebook [https://signalsprocessed.blogspot.de].

The first alternative is scikit.resample. Using it has two problems: it relies on an external (C-language) library called "Secret Rabbit Code" (SRC as in "sample rate conversion", geddit?) aka libsamplerate [http://www.mega-nerd.com/SRC/]. Not a problem in itself, but you need to install it and this step is not automated in pip and friends (as far as I know at time of writing). (Problem for me was that I needed to run it on a platform where I couldn't do this easily.) Problem 2: scikit.resample 0.3.3 is Python 2 only, and the Python 3 version is unofficial (0.4.0-dev, I used [https://github.com/gregorias/samplerate]).

The other one is resampy, which can be found in PyPI or [https://github.com/bmcfee/resampy]. It is easy to install and works with Python 3 (librosa now uses it as preferred resampler over libsamplerate)

TL;DR: if you can install external libs, use scikit.resample (0.4.0-dev). resampy is faster, but not as good, but entirely reasonable. (If MATLAB resample is good enough for you, resampy will serve you just fine and is easier to install)

Now for a comparison with pretty pictures.

Downsampling a sweep

A common task is to convert between the CD sampling rate of 44.1 kHz, and the multiples of 8kHz that originated from the telekon industry. In particular, 48kHz. (Supposedly the rate of 44.1k was chosen in part because is was difficult to convert to 48 kHz; but it is also connected to NTSC timing). So let's have a sweep sampled at 48 kHz. (See the Jupyter Notebook for the Python code)

Sweep 100-23900 Hz, sampled at 48 kHz.

Now convert it to 44.1 kHz using the different libraries, and plot the spectrograms.

Same sweep downsampled to 44.1 kHz using scikit.resample.

scikit.samplerate (or libsamplerate AKA "Secret Rabbit Code") does it well. Almost no aliasing, yet close to Nyquist, and noise levels in line with the input. Note that once the sweep is above Nyquist, the output is basically nil.

The same sweep but now downsampled using scipy.signal.resample.

And this is scipy.signal.resample. Throughout the entire signal there is a high-frequency tone, and there is considerable energy even if the sweep is above Nyquist. This is a result of the frequency-domain processing operating on the entire signal at one. This is NOT suited to audio signals.

The same sweep downsampled using resampy.

Finally, resampy. There are some funny nonlinear effects, but at very low level - for audio applications generally nothing to worry about. Even the aliasing is below -60 dB. So, while not as good as scikit.resample is is useable and certainly better than scipy.signal.resample!

Upsampling an impulse

A way to examine the antialias filter is to upsample an impulse. Here, note a frustrating aspect of all these differing libraries: the interface is completely different beween all of them; and to boot, resampy and scikit.resample give different length outputs for the same upsampling ratio. This means you can't just switch between them using a from foo import resample statement, you'll need a wrapper function if you want to switch freely amongst them (see librosa as an example).

>>> us1 = sk_samplerate.resample(impulse, P/Q, 'sinc_best')
>>> us2 = scipy_signal.resample(impulse, int(len(impulse)*P/Q))
>>> us3 = resampy.resample(impulse, Q, P)
>>> print(us1.shape, us2.shape, us3.shape)
(1087,) (1088,) (1088,)

Frequency response of an impulse upsampled from 44.1 kHz to 48 kHz

'scikit.resample' certainly uses the best antialias filter, with a steep slope, allowing very little energy over Nyquist. 'resampy' is a more gentle filter (lower order, perhaps?) but decent. 'scipy.signal.resample'... well, anyways.

What do the impulse responses look like in time domain?

Time-domain impulse upsampled from 44.1 kHz to 48 kHz

The same section of the impulse, but with samples converted into dB (10*log10 magnitude squared)

Unsurprisingly, 'resampy' is the most compact, which explains the gentler cut-off. 'scikit.resample' is wider, but acceptable. 'scipy.signal.resample' sticks out like a sore thumb.

Speed comparison

>>> %timeit sk_samplerate.resample(sig, Q/P, 'sinc_best')
10 loops, best of 3: 98.1 ms per loop
>>> %timeit scipy_signal.resample(np.float64(sig), int(len(sig)*Q/P))
100 loops, best of 3: 4.52 ms per loop
>>> %timeit resampy.resample(np.float64(sig), P, Q) 
10 loops, best of 3: 40.7 ms per loop

And here is a notable merit of 'scipy.signal.resample'. It is faster by an order of magnitude compared to the other methods. I suspect that if you make sure your signals are of length 2^N, you'll get even faster results, since it'll switch to a FFT instead of a DFT. The other two are probably losing some speed in the passing of data from Python to C - but fundamentally, frequency domain techniques do tend to be fast. So keep that in mind you you need speed.

Conclusion

I'll just repeat the above TL;DR: use scikit.resample (0.4.0-dev) if you can install libsamplerate in a sane place. Else, use resampy, which is faster, but not as good --- but entirely reasonable. Use 'scipy.signal.resample' only if you really need the speed, and be aware of its shortcomings (note the circular assumption!).

FFT-based Overlap-Add FIR filtering in Python

Here is a small Python function I've written (github), that might be useful if you're doing signal processing in Python.  The function implements the classic FFT-based Overlap-Add filtering, potentially saving a heck of a lot of processing time (assuming your filter is of sufficiently high order: usually about 128 taps).

For a thorough explanation of the algorithm, see the Wikipedia article.

I wrote this code as part of a lager ongoing project (gammatone filtering) that I will release eventually.  What I still have to add to this code is the ability to set and save state information (it's simply the part of the response cut off at the end of the function) and the ability to filter complex signals with complex filters (replacing rfft with fft).  Changes made in the latest version.

A simple scikit-learn classifier based on Gaussian Mixture Models (GMM)

When I started switching to Python for my work on CASA, it wasn't entirely clear to me how to use the sklearn GMM (sklearn.mixture.GMM) for classification.  Turned out a bit easier than expected (yay for scikit-learn!), but for others, here is my implementation of a class that behaves like the other classifiers (eg. sklean.svm.SVC).  All you need to decide is how many Gaussians you want to model your data with, and off you go.

Link to Github repo. A Jupyter notebook shows a sample use.

Why isn't something like this in sklearn yet?  Well, turns out someone did propose it already (no surprise) in a much more general way: see this discussion on GitHub.  (I myself was pointed there when I asked about my own code)  My bit of code is far more primitive, but I hope easier to understand.

A bit of Python path hackery

Like many people, I have a directory in $HOME for useful python code I write. Similar to MATLAB's "Documents/MATLAB" directory, I want those files to be easily available if I'm hacking new stuff. The typical way of handling this is to do:

import sys
sys.path.append('/path/to/my/dir')
       
 

but that is a bit ugly, and has the problem that if I'm doing an IPython Notebook (as I often do), this append function gets reevaluated if I reexecute the cell in which I do all my imports (since as I'm edition the notebook and adding stuff, I would do a lot). Besides I now have sys in my namespace.

I could just dump it into ".local/lib/python3.4/site-packages" (too hidden, outside the tree that is synced between machines), or link "Documents/Python" to that path, but then I don't want to have pip install stuff there, and also, there might be times I don't want "Documents/Python" to be in the path.

So here's my solution.  I place a file in ".local/lib/python3.4/site-packages" and ".local/lib/python2.7/site-packages" with a name like "LocalPath.py", and this file contains

       
import sys

myPythonDir = '/home/jthiem/Documents/Python'
if myPythonDir not in sys.path:
    sys.path.append(myPythonDir);
    
 

I can now do

       
import LocalPath
 

and have instant no-fuss access to my homebrew packages.  Now isn't that nice.

The new and improved method of real-time MIDI control of MATLAB (or Python, or ... whatever!)

In an older post on my blog, I had this method to get the state of a MIDI controller into MATLAB (though readable by pretty much anything that can read files).  This method is quite clunky and when I wanted to do something similar again, I re-thought the problem and came up with a better method based on memory mapped files. You can find the code on GitHub.

Basically, I'm creating a file of fixed length (260 bytes, in /tmp), mmap() it and update the contents based on the MIDI stream I'm receiving.  The first 128 bytes are for keys, where each byte is the last seen velocity (0 means "off").  The following 128 bytes are for the cc messages.  Bytes 257 and 258 store the pitch bend, a 14-bit value with MSB in byte 257.  Byte 259 is the last received program change value and 260 is the channel aftertouch.

MATLAB Access

Accessing the values from MATLAB does not require any special functions, it can simply be done using

mm = memmapfile('/tmp/midibroadcast');

then the data can be obtained in real-time from mm.Data.  As described above, mm.Data(1:128) are the last seen key down velocities (0 meaning key is not pressed) and mm.Data(129:256) are controllers 0..127.  The pitch bend value can be obtained as mm.Data(257)*256+mm.Data(258) (note this may be buggy; my cheap controller (KeyRig 25) does not let me set and hold precise pitch bends, so I can't test it properly). mm.Data(259) is the program change and mm.Data(260) shows the current aftertouch amount.

Python Access

From Python the values can be accessed as easily:

import os
import mmap
mfd = os.open('/tmp/midibroadcast', os.O_RDONLY)
mfile = mmap.mmap(mfd, 0, prot=mmap.PROT_READ)

Remembering that Python counts from 0, the controllers can be read in real-time as mfile[128] to mfile[255]: just subtract 1 from the descriptions above.

Any language which can do memory-mapping should be able to do the same, but it should even be possible to read the current state just by re-reading the /tmp/midibroadcast file.

Enjoy!

Using an ARM Chromebook for Scientific (and Academic) Computing

Samsung ARM Chromebook

A couple of months ago, I decided to get myself a Chromebook. The Samsung ARM Chromebook is cheap (to the point of being almost disposable) and it's got an ARM CPU - and for performance at the least possible amount of juice it's hard to beat. I really like the fact that this thing emits no noise that I can detect even in a very quiet room.

But how useful is it, for someone in a standard engineering/academic setting? The answer is that it works well, for me at least - with some special considerations. Especially for the last few weeks, it has been my primary laptop, having been dragged to research cluster meetings and one conference. I will explain the details of a few typical things I do, such as (LaTeX) document editing, intensive numerical computation, etc. Read below the break for details.  

Transplant: yet another bridge between MATLAB and Python, but a good one!

This post is basically just an advertisement for a project done by a Master's student that I'm co-supervising at the moment.  While there are numerous methods already out there to link MATLAB to Python (and rumour has it that the next(?) release will make it easier to call Python from MATLAB), I think Basti's "transplant" (github link) strikes a good balance between simplicity and capability.  Bastian's code is elegant and reliable.  My own contribution has just been a small bug fix, the ability to transfer logical (boolean) matrices, and an attempt to add the ability to capture MATLAB's stdout (Bastian came up with a much better solution).

For me the resulting killer feature is that I can write IPython notebooks that call MATLAB code in a sane way.  Complex code which would take too much effort to convert to Python can be called, then the results can be plotted in the Notebook, which is great when working remotely (a longer post on my workflow is in the works...).  Results become more accessible, with a lot of the complexity hidden away in .m files.

So, check it out and spread the word!

A Sudoku solver in Python

I'm not that great at Python yet.  But practice makes perfect, so just to see if I can do it, I wrote this Sudoku solver in Python (needs v3).  The algorithm is my own though I don't think it is very (or at all) novel, just actually looking for other solvers would spoil the fun!

The implementation is also probably overly byzantine, mostly due to the way the state is stored.  This is done as follows: at the bottom is a list of states for each symbol in a given cell, the states are one of "this symbol is possible here", "this symbol is impossible here", "this symbol is here (set by specification)", and "this symbol is here (set by the solver)".  The last two are really the same for the algorithm but eased my own understanding and debugging - but this adds lines to the code since both states need to be checked in some cases.

The set of states for each cell are then bunched into a list for each row, and then all rows are bunched into a list.  So, to find out if the "4" is possible in column 3, row 7, check if state[3][7][4] is set to 0.

The solver algorithm can now be simply described by three steps: 

1. After setting a cell to some symbol, mark the same symbol as being impossible in cells in the same groups (row, column box). 
2. For all cells, check if there is the case where a single symbol is not impossible.  If so, set the cell to that symbol.
3. For all groups, check where within a group a given symbol is possible in only one cell.  If so, assign that symbol to the open cell.

If in either of step 2 or 3 a cell was set to some symbol, restart from one.  Repeat this until neither of the two last steps sets a new cell to a symbol, at which point the puzzle is solved, it's not completely specified, or there is an ambiguity that I don't handle in this code (I think it is possible for some sort of circular ambiguity to exist).