# HG changeset patch
# User Chris Cannam
# Date 1395746555 0
# Node ID 2b7257e4fc8accc9936ea5d576dfde960f97ab33
# Parent  a91de434feb81d6ec7fffebf3db725402c609d50
use tiledTo, some more notes

diff -r a91de434feb8 -r 2b7257e4fc8a notes/cplcaMT-annotated.m
--- a/notes/cplcaMT-annotated.m	Mon Mar 24 16:31:20 2014 +0000
+++ b/notes/cplcaMT-annotated.m	Tue Mar 25 11:22:35 2014 +0000
@@ -180,7 +180,8 @@
 %%          sum[p,w,f,s] ( P(p,f,s|w,t) P(w,t) )
 %%
 %% (there is also an update equation for x, or P(w|s,p) but we
-%% don't want that as it's the input)
+%% don't want that as it's the input -- one paper proposes an 89th
+%% template to learn the noise component but... not yet)
 
 
 
@@ -202,8 +203,7 @@
 
 	  %% I believe this is equivalent to performing a 553-point
 	  %% FFT of each column of the input (with w{r,k} in the first
-	  %% 545 elements of the first column of that input) and then
-	  %% a 199-point FFT of each row of the result.
+	  %% 545 elements of the first column of that input).
 
 	  %% The output is of course complex.
 
@@ -253,6 +253,9 @@
     
     xbar = x ./ xa;
     xbar = eval( flz);
+
+    %% xbar now contains the result of Eqn 8 in the CMJ paper, Pt(p,f,s|w)
+
     fx = fftn( xbar, wc);
     
     
@@ -283,7 +286,7 @@
                 nh1 = nh1 .* repmat(u{r,k},1,size(h{k},1))';
                 nh = nh + nh1; %% so nh will presumably be 100x5 too
                 
-                nhu = eval( fnh); %% more mystery
+                nhu = eval( fnh); %% more magic
 
 		%% h{k} is 5x100, I'd expect this to be 100x5, I must
 		%% have got something transposed somewhere
diff -r a91de434feb8 -r 2b7257e4fc8a yeti/em.yeti
--- a/yeti/em.yeti	Mon Mar 24 16:31:20 2014 +0000
+++ b/yeti/em.yeti	Tue Mar 25 11:22:35 2014 +0000
@@ -40,9 +40,10 @@
         fold do acc instrument:
             fold do acc note:
                 template = mat.getColumn note data.templates[instrument];
-                w = mat.repeatedHorizontal (mat.width chunk) (mat.newColumnVector template);
-                p = mat.repeatedVertical (mat.height chunk) data.pitches[note];
-                s = mat.repeatedVertical (mat.height chunk) data.sources[instrument][note];
+                resize = mat.tiledTo (mat.size chunk);
+                w = resize (mat.newColumnVector template);
+                p = resize data.pitches[note];
+                s = resize data.sources[instrument][note];
                 mat.sum [acc, mat.entryWiseProduct [w, p, s]];
             done acc [data.ranges[instrument].lowest .. 
                       data.ranges[instrument].highest]
@@ -50,18 +51,20 @@
     mat.entryWiseDivide chunk estimate);
 
 performMaximisation data chunk error =
-   (fold do acc note:
-        fold do acc instrument:
+   (pitches =
+        fold do acc note:
+            fold do acc instrument:
+                // want sum of error * original for all template and instruments
+                // for this pitch, divided by sum of error * original for all
+                // template and instruments for all pitches
+
             template = mat.getColumn note data.templates[instrument];
             w = mat.repeatedHorizontal (mat.width chunk) (mat.newColumnVector template);
             p = mat.repeatedVertical (mat.height chunk) data.pitches[note];
             s = mat.repeatedVertical (mat.height chunk) data.sources[instrument][note];
-
-
-            mat.sum [acc, mat.entryWiseProduct [w, s, error]]
-
+            mat.sum [acc, mat.entryWiseProduct [w, p, s, error]]
         done acc (select do i: inRange data.ranges i note done data.instruments)
-    done (mat.constMatrix epsilon (mat.size chunk)) [data.lowest .. data.highest]);
+    done (mat.constMatrix epsilon  [data.lowest .. data.highest]);
 
 {
     initialise,