view yeti/silvet.yeti @ 26:fbc4011c7693

Now to try to get the matrix version working!

program silvet;

{ prepareTimeFrequency } = load timefreq;
{ loadTemplates, extractRanges } = load templates;

em = load em;
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;

chunkSize = { rows = height, columns = 100 };

emdata = em.initialise ranges templates 88 chunkSize;

eprintln "initialised EM data: overall pitch range \(emdata.lowest) -> \(emdata.highest)";

chunkify cols = 
    if empty? cols then []
    else
       (mat.resizedTo chunkSize
           (mat.fromColumns (take chunkSize.columns cols)))
        :. \(chunkify (drop chunkSize.columns cols));
    fi;

chunks = chunkify columns;

eprintln "we have \(length chunks) chunks of size \(mat.size (head chunks))";

oneIteration emdata chunk n =
   ({ estimate, q } = em.performExpectation emdata chunk;
    newdata = em.performMaximisation emdata chunk q;
    if (n % 6 == 0) then
        \() (plot.plot [ Grid chunk ]);
        \() (plot.plot [ Grid estimate ]);
    fi;
    newdata);

iterations = 6;

var d = emdata;
for [1..iterations] do i:
    eprintln "iteration \(i)...";
    d := oneIteration d (head chunks) i;
done;    

eprintln "done";


/*
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;
*/

();