//
// a beat detector in ChucK
//   and something like a tempo-detector, also in ChucK
// by Tom Lieber
// based on "Tempo and beat analysis of acoustic musical signals"
// by Eric Scheirer
//
// to use, modify the variables in the "configuration" section,
// especially connecting an appropriate sound-producing unit
// generator to "input" and setting a proper value for "length"
//

Gain input;

//
// configuration
//

// connect the sound to beat-detect to "input" and set "length"
SndBuf snd => input;
"clicks.wav" => snd.read;
snd.samples()::samp => dur length;

// number of combs per filter band
50 => int numcombs;

// frequency at which the differentiator samples
200 => float differentiator_freq;

// min size of a differentiator pulse to be a beat
0.03 => float cutoff;

//
// set up input and output unit generators
//

input => Delay inputdelay => dac.chan(0);
0.03::second => inputdelay.max => inputdelay.delay;

Gain output => dac.chan(1);


//
// frequency filterbank
//

// must be 6 (unless you want to recalculate the filters below) (I sure don't)
6 => int bands;

Gain filterbank_outs[bands];

BiQuad bq1[bands];
BiQuad bq2[bands];
BiQuad bq3[bands];

float b1[bands][3];
float b2[bands][3];
float a1[bands][3];
float a2[bands][3];
float bqgain[bands][3];

/*
The coefficients for the BiQuad filters below were generated
by MATLAB R2008b.

This is how I generated each of the filters:
  [Z,P,K] = ellip(6, 3, 40, 0.009070294784580, 'low')
  [Z,P,K] = ellip(3, 3, 40, [0.009070294784580 0.018140589569161])
  [Z,P,K] = ellip(3, 3, 40, [0.018140589569161 0.036281179138322])
  [Z,P,K] = ellip(3, 3, 40, [0.036281179138322 0.07256235827664])
  [Z,P,K] = ellip(3, 3, 40, [0.07256235827664 0.145124716553288])
  [Z,P,K] = ellip(6, 3, 40, 0.145124716553288, 'high')

For each filter, I found the biquad coefficients by doing this:
  [sos,g] = zp2sos(Z,P,K)
  Hd = dfilt.df2tsos(sos,g)
  h = fvtool(Hd)
  click "filter coefficients" in toolbar (far right)
*/

// band 0, low pass, < 200 Hz
-1.994347735657222564498169958824291825294 => b1[0][0];
1.000000000000000444089209850062616169453  => b2[0][0];
-1.988275640269636657109231236972846090794 => a1[0][0];
0.988451492385163499321265589969698339701  => a2[0][0];
0.009935486290798312830618854718522925396  => bqgain[0][0];
-1.998790278411389564894307113718241453171 => b1[0][1];
0.999999999999998889776975374843459576368  => b2[0][1];
-1.995582590576266301951591231045313179493 => a1[0][1];
0.996215113856924294744032977177994325757  => a2[0][1];
1                                          => bqgain[0][1];
-1.999070740216971842073689913377165794373 => b1[0][2];
1.000000000000001110223024625156540423632  => b2[0][2];
-1.99853290274533446080340581829659640789  => a1[0][2];
0.999334599367840192840617419278714805841  => a2[0][2];
1                                          => bqgain[0][2];

// band 1, band pass, 200-400 Hz
 0                                         => b1[1][0];
-1                                         => b2[1][0];
-1.98921720112725353146743145771324634552  => a1[1][0];
 0.990833598456158037848240383027587085962 => a2[1][0];
0.000781585412868364502042561881012261438  => bqgain[1][0];
-1.993039000488599610250162186275701969862 => b1[1][1];
 1.000000000000000222044604925031308084726 => b2[1][1];
-1.991945156231202940233515619183890521526 => a1[1][1];
 0.995016079025449862172081338940188288689 => a2[1][1];
1                                          => bqgain[1][1];
-1.999621583669715096931440712069161236286 => b1[1][2];
 0.999999999999998445687765524780843406916 => b2[1][2];
-1.996512332150246038509067147970199584961 => a1[1][2];
 0.997367570435528372208011660404736176133 => a2[1][2];
1                                          => bqgain[1][2];

// band 2, band pass, 400-800 Hz
 0.000000000000000888178419700125232338905 => b1[2][0];
-1.000000000000000444089209850062616169453 => b2[2][0];
-1.975312391497472397361434559570625424385 => a1[2][0];
 0.981747164593608112603817517083371058106 => a2[2][0];
0.001554408646335501357343966688517866714  => bqgain[2][0];
-1.972241320353174787172179094341117888689 => b1[2][1];
 0.999999999999999333866185224906075745821 => b2[2][1];
-1.977819517440656538198595626454334706068 => a1[2][1];
 0.990063795551808789596748283656779676676 => a2[2][1];
1                                          => bqgain[2][1];
-1.998486948174321931048780243145301938057 => b1[2][2];
 1.000000000000000666133814775093924254179 => b2[2][2];
-1.991324130186075791471012053079903125763 => a1[2][2];
 0.994739921115257774708595661650178954005 => a2[2][2];
1                                          => bqgain[2][2];

// band 3, band pass, 800-1600 Hz
 0.000000000000001554312234475219156593084 => b1[3][0];
-1.000000000000000222044604925031308084726 => b2[3][0];
-1.93831036965515668413218008936382830143  => a1[3][0];
 0.963795601109460720046229198487708345056 => a2[3][0];
0.003091110817045308357953192768263761536  => bqgain[3][0];
-1.890312219058661424853085009090136736631 => b1[3][1];
 1.000000000000001110223024625156540423632 => b2[3][1];
-1.93168574512256663311404736305121332407  => a1[3][1];
 0.980280215995598780409636674448847770691 => a2[3][1];
1                                          => bqgain[3][1];
-1.993957610690330728431263196398504078388 => b1[3][2];
 0.999999999999999444888487687421729788184 => b2[3][2];
-1.975873707417993507817755016731098294258 => a1[3][2];
 0.989490358895407928763177096698200330138 => a2[3][2];
1                                          => bqgain[3][2];

// band 4, band pass, 1600-3200 Hz
-0.000000000000000444089209850062616169453 => b1[4][0];
-1.000000000000000888178419700125232338905 => b2[4][0];
-1.828862749547298438557163535733707249165 => a1[4][0];
 0.92865468802354822486222474253736436367  => a2[4][0];
0.006242925234643366232123451453617235529  => bqgain[4][0];
-1.581699259192773210358495816763024777174 => b1[4][1];
 1.000000000000000666133814775093924254179 => b2[4][1];
-1.771088728416984947600099076225887984037 => a1[4][1];
 0.96137369848423892548794356116559356451  => a2[4][1];
1                                          => bqgain[4][1];
-1.975987663536304284050970636599231511354 => b1[4][2];
 0.9999999999999978905762532122025731951   => b2[4][2];
-1.924943400562922235508267476689070463181 => a1[4][2];
 0.978953907285684987726881445269100368023 => a2[4][2];
1                                          => bqgain[4][2];

// band 5, high pass, > 3200 Hz
-1.969349500292591592653934640111401677132 => b1[5][0];
 1                                         => b2[5][0];
-0.895634416156746659964937862241640686989 => a1[5][0];
 0.483525084917864500511086589540354907513 => a2[5][0];
0.391770128820141338987781409741728566587  => bqgain[5][0];
-1.86055970801238346190586980810621753335  => b1[5][1];
 1.000000000000002220446049250313080847263 => b2[5][1];
-1.679795121826167969558696313470136374235 => a1[5][1];
 0.928625856374035940632438723696395754814 => a2[5][1];
1                                          => bqgain[5][1];
-1.820352888495074594743527995888143777847 => b1[5][2];
 1.000000000000003108624468950438313186169 => b2[5][2];
-1.784021146011410818843501147057395428419 => a1[5][2];
 0.989646570791333113348287042754236608744 => a2[5][2];
1                                          => bqgain[5][2];

for(0 => int i; i < bands; i++) {
    input => bq1[i] => bq2[i] => bq3[i] => filterbank_outs[i];
    
    1 => bq1[i].b0;
    b1[i][0] => bq1[i].b1;
    b2[i][0] => bq1[i].b2;
    1 => bq1[i].a0;
    a1[i][0] => bq1[i].a1;
    a2[i][0] => bq1[i].a2;
    bqgain[i][0] => bq1[i].gain;
    
    1 => bq2[i].b0;
    b1[i][1] => bq2[i].b1;
    b2[i][1] => bq2[i].b2;
    1 => bq2[i].a0;
    a1[i][1] => bq2[i].a1;
    a2[i][1] => bq2[i].a2;
    bqgain[i][1] => bq2[i].gain;
    
    1 => bq3[i].b0;
    b1[i][2] => bq3[i].b1;
    b2[i][2] => bq3[i].b2;
    1 => bq2[i].a0;
    a1[i][2] => bq3[i].a1;
    a2[i][2] => bq3[i].a2;
    bqgain[i][2] => bq3[i].gain;
}

//
// envelope extractors
//

Gain envelope_outs[bands];
for(0 => int i; i < bands; i++) {
    filterbank_outs[i] => FullRect rect => LPF envlpf => HalfRect rect2 => envelope_outs[i];
    envlpf.set( 30, 1 );
}

//
// differentiator + half-wave rectifier
//

Gain differentiator_outs[bands];

fun void differentiator(UGen in, Step out) {
    in => blackhole;
    in.last() => float last;
    while((second / differentiator_freq) => now) {
        in.last() - last => out.next;
        in.last() => last;
    }
}

for(0 => int i; i < bands; i++) {
    Step step => HalfRect rect => differentiator_outs[i];
    0 => step.next; // otherwise, weird click at the beginning
    spork ~ differentiator(envelope_outs[i], step);
}

//
// resonator bank
//

60 => float mintempo;
120 => float maxtempo;

Delay combs[bands][numcombs];
Gain combgains[bands][numcombs];

Gain resonator_outs[bands][numcombs];

fun float map(float val, float frommin, float frommax, float tomin, float tomax) {
    return (val - frommin) / (frommax - frommin) * (tomax - tomin) + tomin;
}

for(0 => int i; i < bands; i++) {
    for(0 => int p; p < numcombs; p++) {
        differentiator_outs[i] => combgains[i][p] => resonator_outs[i][p];
        resonator_outs[i][p] => combs[i][p] => resonator_outs[i][p];
        
        map(p, 0, numcombs-1, mintempo, maxtempo) => float tempo;
        (1.0 / tempo)::minute => dur period;
        period / differentiator_freq => combs[i][p].max => combs[i][p].delay;

        Math.pow(0.5, 1500::ms/period) => float alpha;
        1 - alpha => combgains[i][p].gain;
        alpha => combs[i][p].gain;
    }
}

// only connect the comb filters at the sample rate of the differentiators
fun void resample_resonators() {
    while(second / differentiator_freq => now) {
        for(0 => int band; band < bands; band++)
            for(0 => int comb; comb < numcombs; comb++)
                resonator_outs[band][comb] => blackhole;
        
        1::samp => now;
        
        for(0 => int band; band < bands; band++)
            for(0 => int comb; comb < numcombs; comb++)
                resonator_outs[band][comb] =< blackhole;
    }
}
spork ~ resample_resonators();

fun void print_resonators(int num, int winner, float winval) {
    "" => string s;
    for(0 => int i; i < num; i++) {
        if(i == winner)
            "*" +=> s;
        else
            "-" +=> s;
    }
    
    <<< mintempo, s, maxtempo, "bmp", "(" + winval + ")" >>>;
}

fun void pick_tempo() {
    float values[numcombs];
    int maxcomb;
    
    while(second / differentiator_freq * 10. => now) {
        0 => maxcomb;
        
        for(0 => int comb; comb < numcombs; comb++) {
            0 => values[comb];
            for(0 => int band; band < bands; band++)
                resonator_outs[band][comb].last() +=> values[comb];
            if(values[comb] > values[maxcomb])
                comb => maxcomb;
        }
        
    	print_resonators(numcombs, maxcomb, values[maxcomb]);
    }
}
spork ~ pick_tempo();

//
// beat monitoring
//

Event beat_detected;

100::ms => dur min_beat_interval;

int found_beat[bands];
for(0 => int i; i < bands; i++)
    0 => found_beat[i];

fun void beat_detect(int differentiator) {
    1 => int reset;
    while(second / 200. => now) {
        if(differentiator_outs[differentiator].last() == 0)
            1 => reset;
        
        if(differentiator_outs[differentiator].last() >= cutoff && reset) {
            beat_detected.broadcast();
            0 => reset;
            1 => found_beat[differentiator];
        } else {
            0 => found_beat[differentiator];
        }
    }
}

fun void beat_tick() {
    now - min_beat_interval => time last;
    Impulse i => output;
    while(beat_detected => now) {
        if(now - last >= min_beat_interval) {
            1 => i.next;
            now => last;
        }
    }
}

for(0 => int i; i < bands; i++)
    spork ~ beat_detect(i);
spork ~ beat_tick();

//
// MAUI viewer
//

// If the number of resonators is small, they can be monitored
// live with MAUI sliders by uncommenting the next block of
// code and running this file in miniAudicle.

/* *
MAUI_View control_view;
MAUI_Slider combsliders[numcombs];

fun void goslidersgo(int i) {
    while(200::ms => now)
    {
        0 => float val;
        for(0 => int p; p < bands; p++)
            resonator_outs[p][i].last() +=> val;
        
        val => combsliders[i].value;
        //<<< i, val >>>;
    }
}

for(0 => int i; i < numcombs; i++) {
    "comb" + i => combsliders[i].name;
    combsliders[i].range( 0, 0.1 );
    combsliders[i].size(300, 50);
    combsliders[i].position(0, i * 50);
    control_view.addElement(combsliders[i]);

    spork ~ goslidersgo(i);
}

control_view.display();
/* */

//
// save
//

// Uncomment this next block of code and create a directory
// called "out" to save a copy of the input (dac left channel),
// output (dac right channel), all filter bands, all envelopes,
// and all differentiated signals to wav files in the "out"
// directory with appropriate names.

// This kills performance on my machine but is useful for
// debugging.

/* *
input => WvOut in_out => blackhole;
"out/input.wav" => in_out.wavFilename;

output => WvOut out_out => blackhole;
"out/output.wav" => out_out.wavFilename;

for(0 => int i; i < bands; i++) {
    filterbank_outs[i] => WvOut bank_out => blackhole;
    "out/bank-" + i + ".wav" => bank_out.wavFilename;
    
    envelope_outs[i] => WvOut env_out => blackhole;
    "out/env-" + i + ".wav" => env_out.wavFilename;
    
    differentiator_outs[i] => WvOut diff_out => blackhole;
    "out/diff-" + i + ".wav" => diff_out.wavFilename;
}
/* */

//
// pass time
//

length => now;
<<< "done" >>>;