ASCIIMath creating images

Showing posts with label Matlab. Show all posts
Showing posts with label Matlab. Show all posts

Monday, February 8, 2016

GMM based localizer on custom ASIC model

The model interface hardware with the FPGA in-circuit emulator.
Lukas Gerlach (L) and Christopher Seifert (R) demoing their ASIC model setup, running realtime on a FPGA.

It's always nice to see one's own research code running on real actual hardware with live data rather than just having a simulation in MATLAB.  My colleagues over at the Institut für Mikroelektronische Systeme (IMS) of the Leibnitz Universität in Hannover presented a demo of their hardware at the Hearing4All winter plenary held last week in Soltau.  The code running on the hardware visible is a GMM based localizer originally written by Tobias May, but since heavily modified by myself.  The next step is that we'll write up exactly what we did to make this all work and how well it does - so look out for an article on this in the near future! It's one of the advantages of being at an intengrated cluster; at Hearing4All, pretty much everything related to hearing loss is being investigated: from basic ear physiology, to audiology, models, algorithms, clinical procedures, implants, and new ground-breaking hardware.


Tuesday, February 2, 2016

The Selective Binaural Beamformer: It's out!


After six or so months of going through the peer review gauntlet, our paper on the Selective Binaural Beamformer (or simply SBB) is finally published.  Thanks to all my coauthors (Menno, Daniel, Simon, and Steven) as well as the reviewers (especially reviewer #2, who gave very tough but important feedback) I think this became a very nice paper.  Please go ahead and read it at http://www.asp.eurasipjournals.com/content/2016/1/12 (EURASIP Journal on Advances in Signal Processing, full title "Speech enhancement for multimicrophone binaural hearing aids aiming to preserve the spatial auditory scene"): it's open access, one can read it either at the above address or download a PDF (see the right sidebar on the linked page).  Being open access, it's free and CC-A 4.0 licensed. 

The basic idea behind the algorithm is this: Normally, if using a beamforming algorithm on a binaural hearing aid, the entire auditory image will collapse to the position of the beam direction, that is ALL sound will appear (to the hearing aid user) to originate from the same location.  Various methods have been proposed to fix this - Simon Doclo in particular has done a lot of work on this topic (which is why it was so helpful to have him as coauthor).  My approach to this problem was to take the signal in the STFT domain (ie, the signal is divided into discrete short time frames and narrow frequency bins) and in each "bin" (time-frequency unit) make a decision if the target signal is dominant, or if the background noise is dominant.  In the first case, I use the beamformer output: the signal is enhanced, collapsed, but that's OK - it _should_ be coming from the target direction anyways.  In the second case, I simply use the signal as it comes from the two microphones closest to the ear canal, without processing - hence there is (almost) no difference from the "real" signal reaching the ears.  So, all the benefit of the beamformer without the nasty collapse of the auditory field!

...Well, mostly.  The tricky part is to make a good speech/noise decision (or actually a "target signal"/background noise decision).  But there's a fancy SNR estimator in there, from Adam Kuklasinski (see ref. 19 - I met him in Lisbon where he presented it at EUSIPCO), and that works pretty well.

So if this is the kind of thing that seems interesting to you, read the paper - and I will post some of the sample files (that were used during subjective testing) soonish on my personal homepage.

Friday, November 28, 2014

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!

Sunday, October 5, 2014

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.  

Wednesday, September 24, 2014

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!

Monday, May 5, 2014

Spatial properties of DEMAND

One of the nice plots from the presentation,
which didn't make it into the paper for space
reasons. The plots show the fit of the measured
coherence to the theoretical prediction.
Here is my primary DAGA 2014 paper, where I examine some of the intermicrophone coherence of the DEMAND recordings.  Also, I experiment a little with calibration, using multidimensional scaling.  Not much one can squeeze into two pages.

The presentation was a bit of a bust - they put me in a session more about policy and noise pollution etc. ("Psychoakustik - Lärmschutzpolitik"), so there was not much useful interaction.  Oh well, it happens.

The paper is here, and the presentation slides here.  This might be quite useful if you're using DEMAND.  The question if anyone is actually using DEMAND (other than me) is still open - I would love to hear from anyone who is.

A bit more interesting is a paper written by a M.Sc. student in our lab, on the test results of a hearing aid algorithm we've been working on. The paper is "Erhaltung der räumlichen Wahrnehmung bei Störgeräuschreduktion in Hörgeräten", by Menno Müller, Joachim Thiemann, Daniel Marquardt, Simon Doclo and Steven van de Par, and the method we use to do binaural noise reduction with preservation of spatial awareness is outlined.  A more detailed paper has been submitted to EUSIPCO 2014, and in about 2 or 3 weeks I should find out if that has been accepted.


Friday, July 5, 2013

Multilevel Non-negative Matrix Factorisation à la FASST

The main graphical description of the FASST
NMF procedure from Alexei's paper [1] 

A particularly useful technique I got acquainted with at my old PostDoc position is the use of Non-negative Matrix Factorisation (NMF), in particular how it is being used in the FASST toolbox. While the FASST toolbox is a great comprehensive framework for source separation, its inner workings are a bit obtuse.  At my new group here in Oldenburg, my colleagues are interested in NMF - but rather than ripping it out of FASST, I decided to simply reimplement the key algorithm used (it's not very hard). The final MATLAB code can be found here, with a quick test script here (oh, and you'll need this function for testing as well). Read on for full explanation...

Tuesday, October 16, 2012

DEMAND: Diverse Environments Multichannel Acoustic Noise Database

It's official: our database of recordings with the microphone array described in a previous post is now active and can be found at http://www.irisa.fr/metiss/DEMAND/. Licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License, this is a large chunk of data - so much, we are only putting it on the website in one format, a zip file of 16 mono wav files for each environment.  My preferred format would have been to post it as 16-channel wav files (easier to work with in MATLAB, provided you have loads and loads of RAM) but that decision had to be made.

The original design called for 18 recordings, but one of those was too challenging to do proper around Rennes, and two were unsatisfactory after recording.  So, those will be left for v2.0.  Still, this is a far more diverse set that most other databases I know of - and with more channels.

So if you're into multichannel signal processing... have at it!

Tuesday, November 22, 2011

The MATLAB code for my Thesis

To my shock and surprise, it turns out that someone other than my committee has actually been reading my thesis.  In fact, Mr. Gao Ge (高戈) of Wuhan University has been looking into my work in such detail that a typographical mistake in Eq. 5.2 was discovered: it should read
$f_m = \frac{1000}{4.37}(10^{\frac{ERB_m}{21.4}}-1)$
which he discovered by comparing the thesis with my 2007 INTERSPEECH paper. (For those of you for whom the math is not displaying correctly, that's f_m = 1000/4.37 (10^(ERB_m/21.4)-1): the exponentiation is kinda important...)  I thank Mr. Gao for his vigilance.

As part of his questions, he has asked if he could see the MATLAB code I used to do my simulations.  Now, the code was mostly over a year old, not very clean and distributed all over the file system of my computer, but I decided it's worth the effort to try to collect it all in case someone else wants to continue research in a similar direction.  And I finally sat down and did just that, and as of now, you can find it here or here.  But no, the code is still not very clean.  However, it seems to run - all you need extra AFAIK is the Signal Processing Toolbox.  You also need to supply your own mono, 16 kHz sampled wave file, which you need to point to in line 64 of RunThesisCode.m.

I'm putting this code out there under the Creative Commons Attribution 3.0 Unported License, and make absolutely no claim that this code is useable for any purpose whatsoever, nor that it is an accurate representation of what I used to run the tests.  However, I hope that the code will be useful to understand the ideas I'm trying to explain in my thesis.

Saturday, September 24, 2011

A progress bar for MATLAB scripts

Here is another helpful little MATLAB script written by Prof. Kabal where I made a small modification. The script is basically a wrapper around MATLAB's  waitbar function, displaying a typical progress bar such as pretty much every program in the existence of GUIs has ever had.

Using a single function, you set it up first with
ProgressBar(Frac, Delta, Title)
where Delta or Title are optional: default values are 0.01 and an empty title. Frac is the starting fraction of completion; typically this will be 0.

You update the value with
ProgressBar(Frac)
which will only update the window if the display will change by more than Delta: So, if Delta is 0.05, the display will only be redrawn if it changes the display by more than 5%.

Finally, call the function without arguments to close the pop-out window.

I made a simple modification to check if we really have a screen to display to: if not, the update operation displays a simple percentage on stdout (also controlled by Delta). I made this modification since I tend to run tasks remotely using ssh and screen.  I found the code to check for a display here.

Get the script here: ProgressBar.m

Sunday, September 11, 2011

Real-time control of MATLAB using a MIDI controller

UPDATE 28/11/2014: this post is obsolete, I have created a better method and described it here.

Sometimes it's nice to set variables in MATLAB in real-time, using an actual physical knob.  Like you can find on a MIDI controller.
MIDI controller hooked up to my Mac
MATLAB generally is not geared for real-time data input, and to the best of my knowledge, there is no MIDI library for input.  So, I hacked up a shockingly simple way of getting controller status data into MATLAB: a small helper program passing data through the filesystem.

The helper program does very little: listen to the MIDI bus, if a controller message is seen, write it to a file in the /tmp directory.  In particular, I name the files /tmp/cc/<number>, where <number> is the controller number.  Since controllers can only take values 0–127, the files are one byte long.

The MATLAB side is equally simple.  If you want to read controller value x, open /tmp/cc/x and read the first byte.  Should the read fail, return -1,: let the calling program figure out what to do in that case (eg. read again or use the old value).  The code is simply:

function r = getCCval( n )

fname = sprintf( '/tmp/cc/%d', n );
f = fopen( fname, 'r' );
if f<0
    r = -1;
    return
end

r = fread( f, 1, 'uint8' );
fclose(f);

if isempty(r)
    r = -1;
end

end

The race condition (MATLAB reading the file while it's being written) also returns -1, given the single-byte length of the file that is sufficient.  Again, the calling program should handle it somehow.

The helper program

The actual reading of MIDI messages is messy, platform-dependent code that I'd rather not be writing myself.  So instead I use the rtmidi package, developed just on the other side of the McGill campus; it works on Mac, Linux, and Windows.  The helper program is a simple modification of the program qmidiin.cpp from the tests directory in the rtmidi distribution.

//*****************************************//
//  midicctmp.cpp by Joachim Thiemann 2011
//  modified from
//  qmidiin.cpp by Gary Scavone, 2003-2004.
//
//  read MIDI queue and create files in
//   /tmp/cc with controller state
//
//*****************************************//

#include 
#include 
#include 
#include "RtMidi.h"

#include 
#define SLEEP( milliseconds ) usleep( (unsigned long) (milliseconds * 1000.0) )

#include 

bool done;
static void finish( int ignore ){ done = true; }

void usage( void ) {
  // Error function in case of incorrect command-line
  // argument specifications.
  std::cout << "\nusage: midicctmp \n";
  std::cout << "    where port = the device to use (default = 0).\n\n";
  exit( 0 );
}

int main( int argc, char *argv[] )
{
  RtMidiIn *midiin = 0;
  std::vector message;
  int nBytes;
  double stamp;
  FILE *f;

  // Minimal command-line check.
  if ( argc > 2 ) usage();

  // RtMidiIn constructor
  try {
    midiin = new RtMidiIn();
  }
  catch ( RtError &error ) {
    error.printMessage();
    exit( EXIT_FAILURE );
  }

  // Check available ports vs. specified.
  unsigned int port = 0;
  unsigned int nPorts = midiin->getPortCount();
  if ( argc == 2 ) port = (unsigned int) atoi( argv[1] );
  if ( port >= nPorts ) {
    delete midiin;
    std::cout << "Invalid port specifier!\n";
    usage();
  }

  try {
    midiin->openPort( port );
  }
  catch ( RtError &error ) {
    error.printMessage();
    goto cleanup;
  }

  // Ignore sysex, timing, or active sensing messages.
  midiin->ignoreTypes( true, true, true );

  // Install an interrupt handler function.
  done = false;
  (void) signal(SIGINT, finish);

  // Periodically check input queue.
  std::cout << "Reading MIDI from port ... quit with Ctrl-C.\n";
  while ( !done ) {
    stamp = midiin->getMessage( &message );
    nBytes = message.size();
    if (nBytes==3) {
      char fn[20];
      if ((message[0]>>4)==11) {
        sprintf(fn,"/tmp/cc/%d",(int)message[1]);
        f = fopen(fn,"w");
        fwrite(&message[2],1,1,f);
        fclose(f);
      }
    }

    // If queue empty, sleep for 10 milliseconds.
    if (nBytes==0) SLEEP( 10 );
  }

  // Clean up
 cleanup:
  delete midiin;

  return 0;
}

The code can be compiled the same way the other program in the tests directory are compiled, on my Mac it was g++ -O3 -Wall -I.. -D__MACOSX_CORE__ -o midicctmp midicctmp.cpp Release/RtMidi.o -framework CoreMIDI -framework CoreFoundation -framework CoreAudio.  Note that the program will fail if the directory  /tmp/cc does not exists; just create it manually, or modify the above code to create it at startup.

I have not yet tested it on Linux, but I see no reason for it to not work.  I have an idea for a slightly more fancy event-based method, but that's for another post.

Saturday, February 12, 2011

Playing sounds from MATLAB on Unix

While I honestly don't know if it is still a problem with the latest versions of MATLAB, there has been a problem with the "sound()" function in MATLAB on Macs and Linux platforms with ALSA.  Here is a very simple script that works around the problem by simply creating a temporary wave file, then calling the appropriate command-line function to make the sound.  It's trivial to customize.

In contrast to the real "sound()" function, this one also doesn't block, and returns immediately.  Again, this is trivial to change (by removing the ampersand) but I like this behavior, especially useful when playing longer sound files.

function [ ] = usound( w, fs )
%USOUND Plays sound on ALSA-based linux machines or macs

if nargin < 2
    fs = 16000;
end

filename = ['/tmp/' getenv('USER') '_matplay.wav'];
wavwrite( w, fs, filename );
if ismac
    eval(['!afplay ' filename ' &']);
else
    eval(['!aplay ' filename ' &']);
end

Sunday, February 6, 2011

Voronoi Neighbors

For a current research project, I was considering the problem of finding the neighbors of some given quantization value, basically a method to figure out if two Voronoi regions share a border, where the quantizer is specified simply by the centroids. I wrote a quick and dirty script, but it doesn't work in all cases - the problem is actually quite difficult for vector quantizers.  However, I started without checking the literature, trying to come up with my own approach; now I know better! Below, my crack at the problem, followed by the better, faster and smarter method!

The initial approach
My approach is based on drawing a line between two centroids.  If I find a point on this line that does not get quantized to either of the centroids, I declare the two regions not neighbors.
In this figure, A and B are clearly neighbors, whereas C isn't.  Testing is done using a binary search: I sample the middle of the line, check if that point gets quantized to A or B - if to B, then pick the halfway point from the middle to A, keep going until either a non-A or B point is found or 10 iterations have been done.  Generally for most checks the algorithm stops at the first check (see A and C) only for the neighboring regions will all 10 iterations be computed.

Problem is, this doesn't work very well.  Let's remove D:

Clearly, C now is a neighbor of A, but the straight line is not where the boundary between A and C is located.  So much for that...

The proper way
I found the correct way of solving this problem by looking at how Voronoi diagrams are drawn by MATLAB: using Delaunay triangulation.  Then using MATLAB functions, the lists of neighbors for the centroids can be found by calling delaunay, then TriRep, then the edges method on the resulting structure.  Here is some sample code:


% create a quantizer
V = zeros(2,8);
V(:,1) = [ 0 0 ];
for n=2:8
   V(:,n) = [ sin(2*pi*(n-1)/8) cos(2*pi*(n-1)/8) ];
end
figure;
voronoi(V(1,:),V(2,:));
tri = delaunay(V(1,:),V(2,:));
tr = TriRep(tri,V(1,:)',V(2,:)');


The resulting edges array can now be examined and turned into a list of neighbors for each centroid.  This should scale to higher dimensions, but won't work for 1-D - but there, the neighboring quantizer region problem is trivial anyways.  (Note, I have the Signal Processing Toolbox, so I'm not sure all those functions are in default MATLAB.)

So while my initial quick bit of code didn't work, it did cause me to think about the problem a bit more deeply and I learned a bit about computational geometry.

Thursday, January 13, 2011

Using MATLAB with simple single-threaded scripts on a cluster

McGill University is part of CLUMEQ, and it is relatively easy to gain access to the Krylov computer, a small(ish) system of 300 cores. While finishing up my thesis, I made some use of this cluster to run simulations in MATLAB.  Here are a few notes on efficiently using the system.

Initially, one would assume the best way to utilize the system is to use the Parallel Toolbox, and surely some applications can fit very well to this mold.  In my case, the parallelism within my code was very data-oriented and thus I couldn't get a great speedup.  Instead, I simply used the cluster as a set of easily queued single-thread computers.  A simple script turns a file of MATLAB commands into a series of queue submissions:

cat matcommands | while read line; do
        echo "#!/bin/bash" > temp.sh
        echo "#PBS -l nodes=1:ppn=1,walltime=1:00:00:00" >> temp.sh
        echo "#PBS -V" >> temp.sh
        echo "#PBS -N Bmatlab_s" >> temp.sh
        echo "cd WorkingDirectory" >> temp.sh

        echo "matlab -nojvm -nodisplay -nodesktop -nosplash -r \"" $line '"' >> temp.sh

        msub temp.sh
done

The "matcommands" file is simply a series of single lines of
script(arguments); exit();
script(arguments); exit();
and so on; these will be executed in parallel independently of each other. Very handy if testing a script with varying data.  The script must be written such that the output is saved to a file whose name depends on the arguments.  I then would download all the result files to my desktop and analyze everything there.

This method of using MATLAB on the cluster also works around the restriction that Krylov only has a 16-cpu license for the Parallel toolbox.  Instead, 30 jobs would typically be processed at the same time (which I assume is simply a local policy restriction).