% beat_demo.m
% This script illustrates the use of the beat_* functions to
% implement a simple music audio beat tracker based on dynamic
% programming, as described in:
% D. Ellis 
% <http://www.ee.columbia.edu/~dpwe/pubs/Ellis07-beattrack.pdf "Beat Tracking by Dynamic Programming">
% J. New Music Research 36(1): 51-60, March 2007, 
% DOI: 10.1080/09298210701653344.

% First, load an example sound file
wavfilename = 'mirex06examples/train01.wav';
[d,sr] = wavread(wavfilename);

% Now, calculate the "onset strength" function as the sum of the
% differentiated and half-rectified energy in a Mel-scale
% time-frequency transform:
[onset_fn, osr, sgram, tt, ff] = beat_onset(d,sr);
% osr returns the sampling rate of the frames of onset_fn
% sgram returns the mel spectrogram as an array, with tt and ff
% being the time and frequency labels.
imagesc(tt, ff, sgram); axis xy
plot(tt(1:end-1), onset_fn);
linkaxes([subplot(211), subplot(212)], 'x');

% The dynamic programming approach needs an estimate of global
% tempo, so we calculate one by autocorrelating the onset function,
% applying a bias window, and choosing the biggest peak.
% Optionally, it plots the windowed autocorrelation for us
linkaxes([subplot(211), subplot(212)], 'off');
display = 1;
tempo = beat_tempo(onset_fn, osr, display);

% Now we can run the dynamic programming beat tracker
beats = beat_simple(onset_fn, osr, tempo);

% We have some helper functions: one to plot the beat times on top
% of the Mel spectrogram (which we saved from beat_onset)
beat_plot(beats, '-r', tt, ff, sgram);

% And we can listen to the result; the system-found beats are
% marked by little tone bursts superimposed on the
% original:
beat_play(beats, d, sr);

% We can go straight from soundfile to beats with beat_track, which
% just provides a wrapper around the steps above:
beats = beat_track(wavfilename);

% We can also read in ground truth for the mirex06
% McKinney/Moelants tapping data
tapfilename = 'mirex06examples/train01.txt';
truth = beat_ground_truth(tapfilename);
% By default this only returns the subset of tap records that are
% consistent with the most popular tempo.  Multiple sequences are
% returned in separate cells of a cell array; we can plot and
% listen to them too:
beat_plot(truth{1}, 'xb');
beat_play(truth{1}, d(1:10*sr), sr);
% Can plot all of them, stretched out, too
beat_plot(truth, 'xy', [], [], [], max(ff));

% We can score a beat track against a ground truth by counting how
% many true beats are missed (deletions), and how many
% system-generated beats don't correspond to true beats
% (insertions):
collar = 0.2; % accept a beat within +/-20% of the tempo period
verbose = 1;
[err,ins,del,tru,hh,dd] = beat_score(beats,truth{1},collar,verbose);
% In this case, only the first beat is 'wrong', because the human
% was late to pick up the beat.
% We can score against all the ..
% .. 35 different tap records for this tempo by passing them all to
% the scoring function:
[err,ins,del,tru,hh,dd] = beat_score(beats,truth,collar,verbose);
% We don't line up with all the human taps, but then no single beat
% track ever could because the humans have too much spread.

% We can evaluate the beat tracker against a whole set of examples:
dirname = 'mirex06examples';
files = dir(fullfile(dirname,'*.wav'));
for i = 1:length(files); fnames{i} = fullfile(dirname,files(i).name); end
% Overall average error rate is high, but it's highly variable
% across examples.