Tuesday, July 1, 2008

Sound Card Spectral Analysis

As my first post on this blog I thought I would introduce a very simple function which exploits some of Matlab's high level data acquisition and plotting abilities.

Although this blog's title is 'iheartmatlab' I will also explore areas about Matlab that I don't like so much.

Ok, onto the first code example:

SoundcardSpectralAnalysis

%===============================================================================
% Description : Acquire acoustic data from default system soundcard and plot
% in both time and frequency domain.
%
% Parameters : Fs - Acquisition sample frequency [44000] Hz
% n_bits - Sample size [16] bits
% n_channels - Number of channels to acquire
% from sound card [2]
% update_rate - Polls sound card for data this
% many times per second [5] Hz
%===============================================================================
function soundcardSpectralAnalysis(Fs, n_bits, n_channels, update_rate)
As per the description this function will continuously acquire data from the soundcard at the specified sample frequency, sample size and update rate for however many channels your soundcad supports.

Default Parameters

One thing I don't like about Matlab is the lack of an efficient manner to define the value of default parameters. I have had to resort to the following to initialise my parameters:
% Initialise default parameters if not supplied
if (~exist('Fs', 'var'))
Fs = 44000;
end
if (~exist('n_bits', 'var'))
n_bits = 16;
end
if (~exist('n_channels', 'var'))
n_channels = 2;
end
if (~exist('update_rate', 'var'))
update_rate = 5;
end
The code checks whether a variable exists within the workspace; and if it does not, it creates it with the default value. There are other possible methods that could have been utilised such as:

if (nargin < 4)
update_rate = 5;
if (nargin < 3)
n_channels = 2;
if (nargin < 2)
n_bits = 16;
if (nargin < 1)
sample_frequency = 44000;
end
end
end
end
Although this method uses slightly less lines of code, the level of nesting makes it a bit hard to understand on first glance the purpose of the code. And the approach implemented is at least insensitive to changes in the order of parameters (which is unlikely... but you never know).

My ideal dream solution for default parameters would be similar to C++ (the other language with which I have some experience):

function soundcardSpectralAnalysis(Fs = 44000, n_bits = 16, n_channels = 2, update_rate = 5)

If a parameter does NOT have an assignment against it, then its deemed to be a required parameter.

Initialising Plots

When you are producing plots in Matlab and are going to be continuously updating the contents of the plot, typically overdrawing or updating the previous result (ie updating an line spectrum or vessel track) then I would recommend initializing a plot with your desired visual properties and then updating only the raw data values contained by the plot.

Too often i've seen:
tic
figure
for i = 1:100
cla
plot(i, i, '.')
drawnow
end
toc
Elapsed time is 4.899789 seconds.
A much more efficient alternative:
tic
figure;
point_plot = plot(nan, nan, '.');
for i = 1:100
set(point_plot, 'XData', i);
set(point_plot, 'YData', i);
drawnow
end
toc
Elapsed time is 1.086440 seconds.
So the SoundcardSpectralAnalysis function initialises the time and frequency domain plots and sets up axis bounds / labels:
plot_colors = hsv(n_channels);

% Initialise plots, one above each other in a single figure window
figure;

% Time Domain plot
subplot(2,1,1)
hold on
for i_channel = 1:n_channels
time_domain_plots(i_channel) = plot(nan, nan, ...
'Color', plot_colors(i_channel, :));
end
xlabel('Sample')
ylabel('Counts')

y_max = 2^(n_bits-1);
ylim([-y_max y_max]);

% Frequency Domain plot
subplot(2,1,2)
hold on
for i_channel = 1:n_channels
freq_domain_plots(i_channel) = plot(nan, nan, ...
'Color', plot_colors(i_channel, :));
end
xlabel('Frequency (Hz)')
ylabel('dB re 1 count/sqrt(Hz)')
xlim([0 Fs/2])
ylim([0 70])

Audio Recorder

Now we come to the guts of the data acquisition, the Matlab inbuilt audiorecorder.
% Setup the audiorecorder which will acquire data off default soundcard
audio_recorder = audiorecorder(Fs, n_bits, n_channels);

set(audio_recorder, 'TimerFcn', {@audioRecorderTimerCallback, ...
audio_recorder, ...
time_domain_plots, ...
freq_domain_plots});
set(audio_recorder, 'TimerPeriod', 1/update_rate);
set(audio_recorder, 'BufferLength', 1/update_rate);

% Start the recorder
record(audio_recorder);
The audiorecorder is a nice simple high level abstraction allowing us to retrieve data from a soundcard. However there is limited control as to how often and how much data is retrieved from the sound card.

Ideally for this application we would determine the number of samples to be read from the sound card (sample_frequency/update_rate) and enter a loop reading this many samples each time. However audiorecorder does not contain this functionality. The only option is to specify a timer callback function and its timer period.

In theory this should work perfectly, however in the real world, timers don't get called EXACTLY every 0.2 seconds (or however long their interval). Depending on CPU load and interrupt timings it can vary.

The effect on audiorecoder is that often you will receive too much / too little data when you query it for the currently recorded data. Unfortunately i could not figure out a simple work around. So this function wil occasionally 'skip' and update and include that data into the next frame's analysis. Sorry folks!

Now.. what makes this code work? Simple its this:
set(audio_recorder, 'TimerFcn', {@audioRecorderTimerCallback, ...
audio_recorder, ...
time_domain_plots, ...
freq_domain_plots});
This says that for each time period as configured, execute the 'audioRecorderTimerCallback' function and supply to it these three parameters.

By then starting the recorder you are beginning recording and starting the timer.

Audio Recorder Callback

The callback function we assigned to the audiorecorder will get called approximately ever 1/update_rate seconds. I wont go into too much detail regarding the callback as this post is long enough as is.

The function initially queries the audiorecorder for some of its properties required for accurate frequency analysis.

It then stops the recorder, retrieves the recorded data before starting the recorder again.
% stop the recorder, grab the data, restart the recorder. May miss some data
stop(obj);
data = getaudiodata(obj, data_format);
record(obj);
It then calculates the power spectrum of the audio data and updates the XData and YData fields of the appropriate time / frequency domain plots before forcing a drawing update of the plots.

Error Handling

You might have noticed that the majority of the code in the Audio Recorder callback function is contained within a try block. This is for 2 reasons:

1) If a systematic error is occuring within the function, then it will continue to occur within the function. Each time that function is called. Which is every time the Audio Recorder timer function is triggered. And because at the command line you do not have access to the soundcardSpectralAnalysis workspace, you cannot stop the timer. Your only solution is to close Matlab.
2) You wish to stop the spectral analysis? Well, just closing the figure window will acheive this.

The corresponding catch block does the following:
catch exception
% Stop the recorder and exit
stop(obj)
rethrow(exception)
end
Most important line is "stop(obj)". This stops the recorder and hence the timer function. The error is then rethrown to alert the user to the issue.

In the case of 1) above, you will be able to debug the code if you were extending its functionality withough having to restart Matlab often. In the case of 2), closing the figure window means that when it attempts to update the plot data fields, they do not exist, causing an error to be thrown. Hence the function can be stopped nicely.

Interesting Issue with Audio Recorder Callback Function

The more observant readers might say:
"Why supply "audio_recorder" to the audioRecorderTimerCallback. Particularly as the audio_recorder is the initiater of the callback available as the "obj" variable, and also because audio_recorder is not even used within the callback function whatsoever!"

And I would completely agree with you.

However. If you remove "audio_recorder" from the callback function parameter list, then the timer is never started and never executed:
set(audio_recorder, 'TimerFcn', {@audioRecorderTimerCallback, ...
time_domain_plots, ...
freq_domain_plots});

function audioRecorderTimerCallback(obj, event, ...
time_domain_plots, freq_domain_plots)

!!!!FAILS!!!!
It is particularly puzzling seeing as though audio_recorder is not used anywhere within the callback function. This can be illustrated by calling "clear('audio_recorder')" at the start of the function. It still functions as per expected.

It might be a bug, but possibly me misunderstanding the behaviour of the callback functions. If anyone has any suggestions, feel free to leave them in the comments below.

Anyhow, thanks for reading my long post about a very simple function. Hopefully its highlighted some interesting aspects of Matlab that you may not have been too familiar with.

I aim to provide a new sample each week along with a little explanation.

7 comments:

  1. Rodney -
    First of all, welcome to the MATLAB blogosphere. Your first post looks great, and I'm looking forward to reading what you have to say in subsequent posts. A few thoughts on this post:

    - FWIW, The extra control you would like to see in the audiorecorder (such as specifying behavior based on # samples recorded) is available with the analoginput object, in the Data Acquisition Toolbox.

    - I also find it annoying that there's no way to find/stop a runaway object. I have passed along this enhancement request to our developers.

    - I reproduced the issue you were having with having to pass in the object as a 3rd input argument. I'm not sure exactly what's going on, but it seems related to scoping of the object variable (which actually relates to the above point). Another way I found to avoid this behavior is to pass the audiorecorder object back out to the caller of your main function. Since this keeps the variable alive, the callback fires. Of course, this isn't as appealing a design for your application. I have passed this issue along to our developers.

    Cheers, and good luck with your new blog!

    ReplyDelete
  2. Hey Scott, Thanks for the response!

    - I've used analoginput before hence why I particularly missed its functionality. However i'm trying to do most of my code samples using only the Core Matlab function set to ensure that everyone with a recent copy of Matlab can try the code out and possibly expand it if need be.

    - Your explanation for the audio_recorder scoping actually makes perfect sense. My next code sample gets around this issue by implementing this function as an SoundcardSpectrumAnalyser class. Thus the audiorecorder has ownership and hence it should exist for the duration of the created object.

    ReplyDelete
  3. I like the idea of wrapping this up in a class. This will definitely keep your audiorecorder instance alive.

    Creating a timer callback on an audiorecorder doesn't imply that the timer now owns the instance of that object.

    So, when soundcardSpectralAnalysis goes out of scope, the audiorecorder object (audio_recorder) gets cleaned up, which effectively calls stop(audio_recorder) and destroys the recorder object.

    Really loved the blog post too! I can't wait for the next one.

    -Nick

    ReplyDelete
  4. That was a great post i really loved this. I am trying to use the same idea but in some other way.

    My application uses a simulink model in which i have used the real-time windows target Analog input and analog output blocks. the data from the analog input block goes to the scope where i have saved the incoming data using a structure with time.

    Now what i want to do is that i wana acquire the data and plot the data on my GUI or on a figure just like you have done in real time as it is getting plotted on the scope. But the scope doesn't let me have that data until it is updating the structure with time which i am using to read the data and analyze the data. Can you please help me how can i do it. Or if i need to modify my application then also give me suggestions.
    your help will be really appreciated in this regard.

    Kind Regards Ehtisham

    ReplyDelete
  5. @Thief^Of^Hearts - Sadly I'm not at all familiar with Simulink.

    I had a quick Google and found the following things that may point you in the right direction:

    - A user created Audio Spectrum Analyser model:
    http://www.mathworks.com/matlabcentral/fileexchange/5572

    - This Knowledge Base article suggests that Simulink has a default spectrum analyser block
    http://www.mathworks.com/support/solutions/data/1-2SBZBT.html?product=SL

    ReplyDelete
  6. HI! Thank you for the blog!
    But the code doesn't work in mine,
    When writing

    set(audio_recorder, 'TimerFcn', {@audioRecorderTimerCallback, audio_recorder});

    and later calling

    function audioRecorderTimerCallback(audio_recorder)

    An error occurs saying too many input arguments, have you found this?

    Thank you!
    .......

    ReplyDelete
  7. Hi mr, i have a question concerning the use of this code. i had an audio raw and i want to use this code to calculate the heart rate. how should i apply it? thank u

    ReplyDelete