Chris@12: 
Chris@12: program silvet;
Chris@12: 
Chris@12: { prepareTimeFrequency } = load timefreq;
Chris@13: { loadTemplates, extractRanges } = load templates;
Chris@14: 
Chris@14: em = load em;
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@19: /*
Chris@12: height = if empty? columns then 0 else vec.length (head columns) fi;
Chris@12: 
Chris@13: chunkSize = { rows = height, columns = 100 };
Chris@13: 
Chris@14: emdata = em.initialise ranges templates 88 chunkSize;
Chris@13: 
Chris@14: eprintln "initialised EM data: overall pitch range \(emdata.lowest) -> \(emdata.highest)";
Chris@13: 
Chris@12: chunkify cols = 
Chris@12:     if empty? cols then []
Chris@12:     else
Chris@13:        (mat.resizedTo chunkSize
Chris@13:            (mat.fromColumns (take chunkSize.columns cols)))
Chris@13:         :. \(chunkify (drop chunkSize.columns cols));
Chris@12:     fi;
Chris@12: 
Chris@12: chunks = chunkify columns;
Chris@12: 
Chris@12: eprintln "we have \(length chunks) chunks of size \(mat.size (head chunks))";
Chris@12: 
Chris@14: eprintln "attempting one expectation phase...";
Chris@13: 
Chris@14: error = em.performExpectation emdata (head chunks);
Chris@14: 
Chris@14: eprintln "done, result has dimension \(mat.size error)";
Chris@14: 
Chris@14: eprintln "attempting one maximisation phase...";
Chris@14: 
Chris@14: newP = em.performMaximisation emdata (head chunks) error;
Chris@14: 
Chris@14: eprintln "done";
Chris@14: 
Chris@14: \() (plot.plot [ Grid (head chunks) ]);
Chris@14: \() (plot.plot [ Grid error ]);
Chris@14: 
Chris@14: \() (plot.plot [ Grid newP ]);
Chris@19: */
Chris@19: 
Chris@19: em1data = em1.initialise ranges templates 88;
Chris@19: 
Chris@23: col = head (drop 50 columns); // (drop ((length columns) / 2) 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@20: oneIteration em1data col n =
Chris@19:    (q = em1.performExpectation em1data col;
Chris@20: //    \() (plot.plot [ Caption "Frequency distribution and output of E-step for iteration \(n)", Vector col, Vector q ]);
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@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: normalise v =
Chris@24:    (s = vec.sum v;
Chris@24:     if s > 0 then vec.divideBy s v 
Chris@24:     else v
Chris@24:     fi);
Chris@24: 
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: