[new_beat_detector
ronw@ee.columbia.edu**20061211213548
 
  - added Dan's new and improved beat detector
  - decreased minimum chunk length to 5 frames 
  - ChunkPower became AvgChunkPower
  - TimeCentroid has been disgraced
  - added synth_chunks.pl - script to synthesize each chunk in EDL in a
    separate wave file
  
 
] {
move ./src/com/meapsoft/disgraced/AvgChunkPower.java ./src/com/meapsoft/featextractors/AvgChunkPower.java
move ./src/com/meapsoft/featextractors/ChunkPower.java ./src/com/meapsoft/disgraced/ChunkPower.java
move ./src/com/meapsoft/featextractors/TimeCentroid.java ./src/com/meapsoft/disgraced/TimeCentroid.java
hunk ./Makefile 20
-classpath = $(JARDIR)/getopt-1.0.12.jar:$(JARDIR)/tritonus_misc.jar:$(JARDIR)/tritonus_share.jar:$(JARDIR)/tritonus_mp3.jar:$(JARDIR)/javalayer.jar:$(SRCDIR)
+classpath = $(JARDIR)/getopt-1.0.12.jar:$(JARDIR)/tritonus_misc.jar:$(JARDIR)/tritonus_share.jar:$(SRCDIR)
addfile ./bin/synth_chunks.pl
hunk ./bin/synth_chunks.pl 1
+#!/usr/bin/perl
+
+# Synthesize each chunk in an EDL file in its own wav file.
+#
+# to use do this:
+# synth_chunks.pl input.edl output_base_name
+
+
+$arg_num = 0;
+
+
+foreach $arg (@ARGV)
+{
+	if ($arg_num == 0)
+	{
+		$EDL_input_filename = $arg;
+	}
+	elsif ($arg_num == 1)
+	{
+		$WAV_filename_base = $arg;
+	}
+	$arg_num++;
+}
+
+print "input edl: $EDL_input_filename\n";
+print "output wav: $WAV_filename_base\n";
+
+if (! open INPUT_FILE, $EDL_input_filename)
+{
+	die "can't open file! $!";
+}
+else
+{
+	print "yay, I opened it!\n";
+}
+
+
+$line_number = 0;
+
+while ($line = <INPUT_FILE>)
+{
+
+	$index = index $line, " ";
+	$length = length($line);
+	$end_string = substr($line, $index, $length);
+	$new_string = "0.0$end_string";
+	#print "$new_string";
+
+
+	$EDL_output_filename = "$EDL_input_filename$line_number.edl";
+	$WAV_output_filename = "$WAV_filename_base$line_number.wav";
+
+	if (! open OUTPUT_FILE, ">$EDL_output_filename")
+	{
+		die "can't open output file! $!";
+	}
+
+	print OUTPUT_FILE "$new_string";
+
+	if ($line_number > 0)
+	{
+		#$command = "java -cp MEAPsoft.jar com.meapsoft.Synthesizer -o $WAV_output_filename $EDL_output_filename\n";
+		system "java -cp MEAPsoft.jar com.meapsoft.Synthesizer -o $WAV_output_filename $EDL_output_filename";
+		print $command;
+	}
+
+	unlink($EDL_output_filename);
+
+	$line_number++;
+}
+
hunk ./src/com/meapsoft/DSP.java 996
+
+    /**
+     * Take the natural exp of each element of a.
+     */
+    public static double[] exp(double[] a)
+    {
+        double[] y = new double[a.length];
+
+        for(int x = 0; x < y.length; x++)
+            y[x] = Math.exp(a[x]);
+
+        return y;
+    }
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 33
- * 3 main problems:
+ * 3 stages:
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 39
- * 2. Have a way to choose a local period
- *    Here, we calculate running autocorrelations and choose a good 
- *    lag based on the YIN trick to discount the zero-order lag, 
- *    and some heuristics to choose among remaining peaks
- *     (which may not be that good).
- *    We also tried taking the short-time Fourier transform
- *    It's not clear how best to enforce local continuity - 
- *    want dynamic programming really.
- *
- * 3. Choose the best phase alignment (downbeat) given a period and 
- *    an approximate location.  Here, we cross-correlate a template 
- *    based on the current period with the envelope to find the 
- *    best match (within +/- 1/2 period, over a -2.5:2.5 beat window).
- *    We also tried taking the maxima of the best STFT component
+ * 2. Estimate the global tempo
+ *    We autocorrelate the entire onset envelope and choose the best 
+ *    peak in the acceptable range (around 100-250 bpm).
+ * 
+ * 3. Go through the onset envelope choosing the best set of times 
+ *    that both lie near a lot of maxima in the envelope and are 
+ *    spaced apart by the global tempo period.  We do this with 
+ *    dynamic programming, choosing a best predecessor for every beat,
+ *    then finally tracing back the beat that has the best total 
+ *    score at the end.
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 64
+    // DP search constants
+    private double tightness = 6.0;	// strictness to tempo 
+    private double alpha = 0.9;		// how much history is weighted vs. local events
+
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 81
-    // The final algo:
-    //   - choose an overall period based on the global autocorrelation
-    //     and stick to it throughout
-    //   - shift the local beats based on best cross-correlation
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 120
+	int pd = 0;
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 136
-        //if length(startpd) == 0
-        //   // I think this means there were no local max in the autoco, unlikely
-        //   error('no global period found');
-        //System.out.println(sgsrate);
-        
-        //System.out.println("Global bpm = "+(60/(startpd/sgsrate))
-                           //+"   startpd = "+startpd);
+	pd = startpd;
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 138
-        // guess first beat as best one around
-        int startcol = (int)Math.round((.5+stabwin)*startpd);
-        int nextbeat = DSP.argmax(DSP.slice(mm, startcol+1, startcol+startpd));
-        nextbeat = nextbeat + startcol;
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 139
-        int pd = startpd;
-        int[] mmlm = DSP.find(localmax(mm));
+	//% Smooth beat events
+	//templt = exp(-0.5*(([-pd:pd]/(pd/32)).^2));
+	double[] templt = DSP.exp(DSP.times(DSP.power(DSP.times(DSP.range(-pd,pd),32.0/pd),2.0),-0.5));
+	//localscore = conv(templt,onsetenv);
+	//localscore = localscore(round(length(templt)/2)+[1:length(onsetenv)]);
+	double[] localscore = DSP.conv(templt,fmm);
+	localscore = DSP.slice(localscore, templt.length/2, templt.length/2+fmm.length-1);
+        //DSP.imagesc(localscore, "localscore");
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 148
-        // cross-correlate with template based on current period
-        // use e.g. +/- stabwin beats
-        int tpd = (int)Math.round((.5+stabwin)*pd);
-        double[] ii = DSP.range(-tpd,tpd);
-        double[] templt = DSP.times(
-            DSP.power(DSP.times(DSP.plus(DSP.cos(DSP.times(ii, 2*Math.PI/pd)), 1), 0.5), 8), 
-            DSP.hanning(ii.length));
-        int hpd = (int)Math.round(pd/2);
+	//% DP version:
+	//% backlink(time) is index of best preceding time for this point
+	//% cumscore(time) is total cumulated score to this point
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 152
-        //double[][] fmmxs = new double[100][Math.round(3*pd/2)];
-        //double[][] xcrs = new double[100][Math.round(pd/2)];
-        //int cb = 0;
-        while(nextbeat + (.5+stabwin)*pd < mm.length)
-        {
-            double[] fmmx = DSP.slice(fmm, nextbeat-tpd, nextbeat+tpd);
-            double[] xcr = DSP.xcorr(fmmx, DSP.times(templt,DSP.max(fmmx)), hpd);
-            int mmx = DSP.argmax(DSP.times(DSP.todouble(localmax(xcr)),xcr));
-            if(mmx == -1)
-            {
-                // xcorr was concave! Use midpoint
-                //System.out.println("t = "+ nextbeat/sgsrate +" sec: xcorr was concave! Use midpoint.");
-                mmx = hpd+1;
-            } 
-            mmx = mmx + nextbeat - hpd - 1;
+	//backlink = zeros(1,length(localscore));
+	//cumscore = zeros(1,length(localscore));
+
+	int[] backlink = new int[localscore.length];
+	double[] cumscore = new double[localscore.length];
+
+	//% search range for previous beat
+	//prange = round(-2*pd):-round(pd/2);
+	int prangemin = -2*pd;
+	int prangemax = -pd/2;
+
+	//% Skewed window
+	//txwt = exp(-0.5*((tightness*log(prange/-pd)).^2));
+	double[] txwt = DSP.exp(DSP.times(DSP.power(DSP.times(DSP.log(DSP.times(DSP.range(prangemin,prangemax),-1.0/pd)),tightness),2.0), -0.5));
+
+	//starting = 1;
+	int starting = 1;
+	double maxlocalscore = DSP.max(localscore);
+
+	//for i = 1:length(localscore)
+	for(int i = 0; i < localscore.length; ++i) {
+  
+	    //timerange = i + prange;
+	    double[] scorecands = new double[txwt.length];
+	    double[] valvals;
+
+	    //% Are we reaching back before time zero?
+	    //zpad = max(0, min(1-timerange(1),length(prange)));
+	    if(i+prangemin < 0) {
+		int valpts = 0;
+		for(int j = 0; j < scorecands.length; ++j) scorecands[j]=0;
+		if(i+prangemax >=0) {
+		    valpts = i + prangemax + 1;
+		    valvals = DSP.times(DSP.slice(txwt,txwt.length-valpts,txwt.length-1),DSP.slice(cumscore,i+prangemax-valpts+1,i+prangemax));
+		    for(int j = 0; j < valpts; ++j) scorecands[scorecands.length-valpts+j] = valvals[j];
+		}
+	    } else {
+	    //% Search over all possible predecessors and apply transition 
+	    //% weighting
+	    //scorecands = txwt .* [zeros(1,zpad),cumscore(timerange(zpad+1:end))];
+		scorecands = DSP.times(txwt,DSP.slice(cumscore,i+prangemin,i+prangemax));
+	    }
+
+	    //% Find best predecessor beat
+	    //[vv,xx] = max(scorecands);
+	    double vv = DSP.max(scorecands);
+	    int xx = DSP.argmax(scorecands);
+	    //% Add on local score
+	    //cumscore(i) = alpha*vv + (1-alpha)*localscore(i);
+	    cumscore[i] = alpha*vv + (1-alpha)*localscore[i];
+
+	    // % special case to catch first onset
+	    //if starting == 1 & localscore(i) < 0.01*max(localscore);
+	    if ( starting == 1 && localscore[i] < 0.01*maxlocalscore) {
+		//backlink(i) = -1;
+		backlink[i] = -1;
+	    } else {
+		//backlink(i) = timerange(xx);
+		backlink[i] = i+prangemin+xx;
+		//% prevent it from resetting, even through a stretch of silence
+		starting = 0;
+	    }
+	}
+
+	//%%%% Backtrace
+
+	//% Cumulated score is stabilized to lie in constant range, 
+	//% so just look for one near the end that has a reasonable score
+	//medscore = median(cumscore(localmax(cumscore)));
+	double medscore = DSP.median(DSP.subsref(cumscore, localmax(cumscore)));
+	//bestendx = max(find(cumscore .* localmax(cumscore) > 0.5*medscore));
+	int bestendx = 0;
+	int jj = cumscore.length-2;
+	while(jj > 0 && bestendx == 0) {
+	    if (cumscore[jj] > cumscore[jj-1] && cumscore[jj] >= cumscore[jj+1] && cumscore[jj] > 0.5*medscore) {
+		bestendx = jj;
+	    }
+	    --jj;
+	}
+
+	//b = bestendx;
+
+	int nbeats = 0;
+	int[] tmplinks = new int[cumscore.length];
+	tmplinks[0] = bestendx;
+	//while backlink(b(end)) > 0
+	int bb;
+	while ((bb = backlink[tmplinks[nbeats]]) > 0) {
+	    //b = [b,backlink(b(end))];
+	    ++nbeats;
+	    tmplinks[nbeats] = bb;
+	}
+
+	//b = fliplr(b);
+	++nbeats;
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 248
-            // but force it to be a local max of mm, if one nearby
-            double[] tmp = DSP.abs(DSP.minus(mmx, DSP.todouble(mmlm)));
-            int ix = DSP.argmin(tmp);
-            //if(DSP.min(tmp) < pd/4)
-            if(DSP.min(tmp) < pd/8)
-                mmx = mmlm[ix];
-            //else
-                //System.out.println("t = "+nextbeat/sgsrate+" sec: no good local max at "+mmx);
+	//% return beat times in secs
+	//b = b / sgsrate;
hunk ./src/com/meapsoft/DpweBeatOnsetDetector.java 251
-            int thisbeat = mmx;
-            
-            nextbeat = thisbeat + pd;
-            
-            //System.out.println("thisbeat = "+thisbeat+" nextbeat = "+nextbeat+" pd = "+pd);
+	for (int j = 0; j < nbeats; ++j) {
+	    int thisbeat = tmplinks[nbeats-1-j];
+	    //System.out.println("thisbeat = "+thisbeat+" nextbeat = "+nextbeat+" pd = "+pd);
hunk ./src/com/meapsoft/DpweOnsetDetector.java 46
-    protected double silenceThresh = -80; 
+    protected double silenceThresh = -40; 
hunk ./src/com/meapsoft/DpweOnsetDetector.java 53
+
hunk ./src/com/meapsoft/DpweOnsetDetector.java 60
-    private double[] boost;
+    // dpwe debug 2006-12-11
+    //static final int nDebugFrames = 1500;
+    //private double[][] melFrames = new double[nDebugFrames][nmel];
+    //private int nMelFrames = 0;
+
+    //private double[] boost;
hunk ./src/com/meapsoft/DpweOnsetDetector.java 80
-        double[] freqs = DSP.times(DSP.range(0,swin/2), sr/swin);
-        double[] lfboost = DSP.times(DSP.log10(DSP.max(DSP.plus(DSP.rdivide(DSP.minus(200, freqs), 200.0), 1), 1)), 10.0);
-        // high frequency boost for noisy drums
-        //double[] hfboost = DSP.times(DSP.log10(DSP.max(DSP.plus(DSP.rdivide(DSP.minus(freqs, 6000), 6000), 1), 1)), 10.0);
-        double[] hfboost = new double[lfboost.length];
-        for(int x = 0; x < hfboost.length; x++)
-        {
-            if(freqs[x] < 6000)
-                hfboost[x] = 0;
-            else
-                hfboost[x] = 3;
-        }
-        boost = DSP.plus(lfboost, hfboost);
+        //double[] freqs = DSP.times(DSP.range(0,swin/2), sr/swin);
+	//        double[] lfboost = DSP.times(DSP.log10(DSP.max(DSP.plus(DSP.rdivide(DSP.minus(200, freqs), 200.0), 1), 1)), 10.0);
+        // // high frequency boost for noisy drums
+        // //double[] hfboost = DSP.times(DSP.log10(DSP.max(DSP.plus(DSP.rdivide(DSP.minus(freqs, 6000), 6000), 1), 1)), 10.0);
+        // double[] hfboost = new double[lfboost.length];
+        //for(int x = 0; x < hfboost.length; x++)
+        //{
+        //    if(freqs[x] < 6000)
+        //        hfboost[x] = 0;
+        //    else
+        //        hfboost[x] = 3;
+        //}
+        //boost = DSP.plus(lfboost, hfboost);
hunk ./src/com/meapsoft/DpweOnsetDetector.java 113
-        if(newestFrame <= numFrames)
+        if(newestFrame <= numFrames && newestFrame != -1)
hunk ./src/com/meapsoft/DpweOnsetDetector.java 116
-            double[] D = stft.getFrame(newestFrame);
-            stft.setFrame(newestFrame, DSP.plus(D, boost)); 
+            //double[] D = stft.getFrame(newestFrame);
+            //stft.setFrame(newestFrame, DSP.plus(D, boost)); 
hunk ./src/com/meapsoft/DpweOnsetDetector.java 126
+
+	    // dpwe debug 2006-12-11
+// 	    if (nMelFrames < nDebugFrames) {
+// 		for (int j = 0; j < nmel; ++j) {
+// 		    melFrames[nMelFrames][j] = melFrame[j];
+// 		}
+// 		++nMelFrames;
+// 		if (nMelFrames == nDebugFrames) {
+// 		    DSP.imagesc(melFrames, "melFrames");
+// 		}
+// 	    }
hunk ./src/com/meapsoft/DpweOnsetDetector.java 150
-
+        }
+        else
+        {
hunk ./src/com/meapsoft/DpweOnsetDetector.java 155
-            if(newestFrame == numFrames)
-            {
-                checkOnsets();
-            }
+            checkOnsets();
hunk ./src/com/meapsoft/MEAPUtil.java 439
+        AudioFormat baseFormat = stream.getFormat();
+        AudioFormat decodedFormat = 
+            new AudioFormat(AudioFormat.Encoding.PCM_SIGNED,
+                            baseFormat.getSampleRate(),
+                            16,
+                            baseFormat.getChannels(),
+                            baseFormat.getChannels() * 2,
+                            baseFormat.getSampleRate(),
+                            bigEndian);
+        stream = AudioSystem.getAudioInputStream(decodedFormat, stream);
+        //System.out.println(stream.getFormat());
+
hunk ./src/com/meapsoft/MEAPUtil.java 464
-        //System.out.println(stream.getFormat());
+        //System.out.println("final: "+stream.getFormat());
hunk ./src/com/meapsoft/STFT.java 56
-  public float samplingRate = 44100;  
+  public float samplingRate = MEAPUtil.samplingRate;  
hunk ./src/com/meapsoft/STFT.java 88
-  public void start() {
+  // returns the total number of bytes read
+  public long start() {
hunk ./src/com/meapsoft/STFT.java 98
+      int totalBytesRead = 0;
hunk ./src/com/meapsoft/STFT.java 110
+              totalBytesRead += bytesRead;
hunk ./src/com/meapsoft/STFT.java 113
-              return;
+              return totalBytesRead;
hunk ./src/com/meapsoft/STFT.java 154
-          
+
hunk ./src/com/meapsoft/STFT.java 159
+
+      // let the frame listeners know that we're done reading:
+      notifyListeners(-1);
+      
+      return totalBytesRead;
hunk ./src/com/meapsoft/SegmentExtractor.java 43
-    OnsetDetector onsDetect;
hunk ./src/com/meapsoft/SegmentExtractor.java 45
+
hunk ./src/com/meapsoft/SegmentExtractor.java 58
-    public SegmentExtractor(STFT stft, OnsetDetector od, String sfn, FeatFile of, BoundedRangeModel brm)
+    public SegmentExtractor(STFT stft, String sfn, FeatFile of, BoundedRangeModel brm)
hunk ./src/com/meapsoft/SegmentExtractor.java 61
-        onsDetect = od;
hunk ./src/com/meapsoft/SegmentExtractor.java 74
+
hunk ./src/com/meapsoft/Segmenter.java 221
-        extractor = new SegmentExtractor(stft, detector, fn, outFile, progress);
+        extractor = new SegmentExtractor(stft, fn, outFile, progress);
hunk ./src/com/meapsoft/disgraced/ChunkPower.java 23
-package com.meapsoft.featextractors;
+package com.meapsoft.disgraced;
hunk ./src/com/meapsoft/disgraced/ChunkPower.java 26
+import com.meapsoft.featextractors.FeatureExtractor;
hunk ./src/com/meapsoft/featextractors/AvgChunkPower.java 23
-package com.meapsoft.disgraced;
+package com.meapsoft.featextractors;
hunk ./src/com/meapsoft/featextractors/AvgChunkPower.java 26
-import com.meapsoft.featextractors.FeatureExtractor;
hunk ./src/com/meapsoft/featextractors/AvgChunkPower.java 50
-		return "I am a generic FeatureExtractor";
+		return "Computes the average power in each chunk.";
}
