diff yeti/silvet_onecolumn.yeti @ 27:9ec18d453889

Split single-column and matrix versions

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/yeti/silvet_onecolumn.yeti	Mon Mar 31 17:06:24 2014 +0100
@@ -0,0 +1,104 @@
+
+program silvet_onecolumn;
+
+{ prepareTimeFrequency } = load timefreq;
+{ loadTemplates, extractRanges } = load templates;
+
+em1 = load em_onecolumn;
+
+mat = load may.matrix;
+vec = load may.vector;
+plot = load may.plot;
+
+templates = loadTemplates ();
+
+ranges = extractRanges templates;
+
+eprintln "\nWe have \(length (keys templates)) instruments:";
+for (sort (keys templates)) do k:
+    eprintln " * \(k) \(mat.size templates[k]) range \(ranges[k].lowest) -> \(ranges[k].highest)";
+done;
+eprintln "";
+
+columns = prepareTimeFrequency "test.wav";
+
+height = if empty? columns then 0 else vec.length (head columns) fi;
+
+em1data = em1.initialise ranges templates 88;
+
+col = head (drop 50 columns); // (drop ((length columns) / 2) columns);
+
+\() (plot.plot [ Caption "Source frequency distribution", Vector col ]);
+
+\() (plot.plot [ Caption "Source distribution beforehand", Grid em1data.sources]);
+
+normalise v =
+   (s = vec.sum v;
+    if s > 0 then vec.divideBy s v 
+    else v
+    fi);
+
+oneIteration em1data col n =
+   ({ estimate, q } = em1.performExpectation em1data col;
+    newdata = em1.performMaximisation em1data col q;
+    if (n % 6 == 0) then
+        \() (plot.plot [ Caption "Pitch distribution before and after M-step update for iteration \(n)", Vector (em1data.pitches), Vector (newdata.pitches) ]);
+        \() (plot.plot [ Caption "Source distribution after M-step update for iteration \(n)", Grid newdata.sources ]);
+        \() (plot.plot [ Caption "Q function for E-step iteration \(n)", Vector q ]);
+        \() (plot.plot [ Caption "Estimate from E-step iteration \(n) against original source distribution, and difference between them", Vector estimate, Vector (normalise col), Vector (vec.subtract estimate (normalise col)) ]);
+    fi;
+    newdata);
+
+iterations = 12;
+
+var d = em1data;
+
+for [1..iterations] do i:
+    d := oneIteration d col i;
+done;    
+
+var sounding = [];
+
+println "pitch distribution: \(vec.list d.pitches)";
+
+for [d.lowest .. d.highest] do p:
+    if (vec.at d.pitches p) > 0.05 then
+        sounding := sounding ++ [p];
+    fi;
+done;
+
+println "Sounding: \(sounding)";
+
+toNote p = (array ["A","A#","B","C","C#","D","D#","E","F","F#","G","G#"])[p % 12];
+
+println "Notes: \(map toNote sounding)";
+
+var instruments = [];
+for sounding do p:
+    var best = 0;
+    var bestp = 0;
+    for [0..d.instCount-1] do i:
+        if mat.at d.sources p i > bestp then
+            bestp := mat.at d.sources p i;
+            best := i;
+        fi;
+    done;
+    if bestp > 0 then
+        instruments := instruments ++ [best];
+    else
+        instruments := instruments ++ [-1];
+    fi;
+done;
+
+println "Instruments: \(map do i: (d.instruments[i]) done instruments)";
+
+if not (empty? sounding) then
+   p = head sounding;
+   i = head instruments;
+   w = mat.getColumn p d.templates[i];
+   \() (plot.plot [ Vector w, Vector (normalise col), Caption "Template for instrument \(d.instruments[i]), pitch \(p), against normalised source distribution" ]);
+fi;
+
+();
+
+