Chris@12: 
Chris@27: program silvet_onecolumn;
Chris@12: 
Chris@12: { prepareTimeFrequency } = load timefreq;
Chris@13: { loadTemplates, extractRanges } = load templates;
Chris@14: 
Chris@19: em1 = load em_onecolumn;
Chris@12: 
Chris@12: mat = load may.matrix;
Chris@12: vec = load may.vector;
Chris@14: plot = load may.plot;
Chris@12: 
Chris@12: templates = loadTemplates ();
Chris@12: 
Chris@13: ranges = extractRanges templates;
Chris@13: 
Chris@13: eprintln "\nWe have \(length (keys templates)) instruments:";
Chris@13: for (sort (keys templates)) do k:
Chris@14:     eprintln " * \(k) \(mat.size templates[k]) range \(ranges[k].lowest) -> \(ranges[k].highest)";
Chris@13: done;
Chris@12: eprintln "";
Chris@12: 
Chris@12: columns = prepareTimeFrequency "test.wav";
Chris@12: 
Chris@12: height = if empty? columns then 0 else vec.length (head columns) fi;
Chris@12: 
Chris@19: em1data = em1.initialise ranges templates 88;
Chris@19: 
Chris@28: col = head (drop 50 columns);
Chris@19: 
Chris@20: \() (plot.plot [ Caption "Source frequency distribution", Vector col ]);
Chris@19: 
Chris@22: \() (plot.plot [ Caption "Source distribution beforehand", Grid em1data.sources]);
Chris@20: 
Chris@25: normalise v =
Chris@25:    (s = vec.sum v;
Chris@25:     if s > 0 then vec.divideBy s v 
Chris@25:     else v
Chris@25:     fi);
Chris@25: 
Chris@20: oneIteration em1data col n =
Chris@25:    ({ estimate, q } = em1.performExpectation em1data col;
Chris@19:     newdata = em1.performMaximisation em1data col q;
Chris@21:     if (n % 6 == 0) then
Chris@20:         \() (plot.plot [ Caption "Pitch distribution before and after M-step update for iteration \(n)", Vector (em1data.pitches), Vector (newdata.pitches) ]);
Chris@22:         \() (plot.plot [ Caption "Source distribution after M-step update for iteration \(n)", Grid newdata.sources ]);
Chris@25:         \() (plot.plot [ Caption "Q function for E-step iteration \(n)", Vector q ]);
Chris@25:         \() (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)) ]);
Chris@20:     fi;
Chris@19:     newdata);
Chris@19: 
Chris@21: iterations = 12;
Chris@19: 
Chris@20: var d = em1data;
Chris@20: 
Chris@20: for [1..iterations] do i:
Chris@20:     d := oneIteration d col i;
Chris@20: done;    
Chris@14: 
Chris@21: var sounding = [];
Chris@21: 
Chris@22: println "pitch distribution: \(vec.list d.pitches)";
Chris@21: 
Chris@21: for [d.lowest .. d.highest] do p:
Chris@21:     if (vec.at d.pitches p) > 0.05 then
Chris@21:         sounding := sounding ++ [p];
Chris@21:     fi;
Chris@21: done;
Chris@21: 
Chris@21: println "Sounding: \(sounding)";
Chris@21: 
Chris@21: toNote p = (array ["A","A#","B","C","C#","D","D#","E","F","F#","G","G#"])[p % 12];
Chris@21: 
Chris@21: println "Notes: \(map toNote sounding)";
Chris@21: 
Chris@21: var instruments = [];
Chris@21: for sounding do p:
Chris@22:     var best = 0;
Chris@21:     var bestp = 0;
Chris@22:     for [0..d.instCount-1] do i:
Chris@22:         if mat.at d.sources p i > bestp then
Chris@22:             bestp := mat.at d.sources p i;
Chris@21:             best := i;
Chris@21:         fi;
Chris@21:     done;
Chris@21:     if bestp > 0 then
Chris@21:         instruments := instruments ++ [best];
Chris@21:     else
Chris@22:         instruments := instruments ++ [-1];
Chris@21:     fi;
Chris@21: done;
Chris@21: 
Chris@22: println "Instruments: \(map do i: (d.instruments[i]) done instruments)";
Chris@21: 
Chris@24: if not (empty? sounding) then
Chris@24:    p = head sounding;
Chris@24:    i = head instruments;
Chris@24:    w = mat.getColumn p d.templates[i];
Chris@24:    \() (plot.plot [ Vector w, Vector (normalise col), Caption "Template for instrument \(d.instruments[i]), pitch \(p), against normalised source distribution" ]);
Chris@24: fi;
Chris@21: 
Chris@14: ();
Chris@14: 
Chris@14: