Chris@16: function [w,h,z,xa] = hnmf( x, K, R, iter, sh, sz, w, h, pl, lh, lz)
Chris@16: % function [w,h,z,xa] = hnmf( x, K, R, iter, sh, sz, w, h, pl, lh, lz)
Chris@16: %
Chris@16: % Perform Multi-source NMF
Chris@16: %
Chris@16: % Inputs:
Chris@16: %  x     input distribution
Chris@16: %  K     number of components
Chris@16: %  R     number of sources
Chris@16: %  iter  number of EM iterations [default = 100]
Chris@16: %  sh    sparsity of h
Chris@16: %  sz    sparsity of z
Chris@16: %  w     initial value of w
Chris@16: %  h     initial value of h
Chris@16: %  pl    plot flag
Chris@16: %  lh    update h flag
Chris@16: %  lz    update z flag
Chris@16: %
Chris@16: % Outputs:
Chris@16: %  w   spectral bases
Chris@16: %  h   component activation
Chris@16: %  z   source activation per component
Chris@16: %  xa  approximation of input
Chris@16: %
Chris@16: % Emmanouil Benetos 2011
Chris@16: 
Chris@16: 
Chris@16: % Get sizes
Chris@16: [M,N] = size( x);
Chris@16: sumx = sum(x);
Chris@16: 
Chris@16: % Default training iterations
Chris@16: if ~exist( 'iter')
Chris@16:     iter = 100;
Chris@16: end
Chris@16: 
Chris@16: % Default plot flag
Chris@16: if ~exist( 'pl')
Chris@16:     pl = 1;
Chris@16: end
Chris@16: 
Chris@16: % Initialize
Chris@16: if ~exist( 'w') || isempty( w)
Chris@16:     w = rand( M, R, K);
Chris@16: end
Chris@16: for r=1:R
Chris@16:     for k=1:K
Chris@16:         w(:,k,r) = w(:,k,r) ./ sum(w(:,k,r));
Chris@16:     end;
Chris@16: end;
Chris@16: if ~exist( 'h') || isempty( h)
Chris@16:     h = rand( K, N);
Chris@16: end
Chris@16: n=1:N;
Chris@16: h(:,n) = repmat(sumx(n),K,1) .* (h(:,n) ./ repmat( sum( h(:,n), 1), K, 1));
Chris@16: if ~exist( 'z') || isempty( z)
Chris@16:     z = rand( R, K, N);
Chris@16: end
Chris@16: for k=1:K
Chris@16:     for n=1:N
Chris@16:         z(:,k,n) = z(:,k,n) ./ sum(z(:,k,n));
Chris@16:     end;
Chris@16: end;
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: % Iterate
Chris@16: for it = 1:iter
Chris@16:     
Chris@16:     % E-step
Chris@19:     zh = z .* permute(repmat(h,[1 1 R]),[3 1 2]); %% z is the source activation distribution, h the component (pitch) activation
Chris@16:     xa=eps;
Chris@16:     for r=1:R
Chris@16:         for k=1:K
Chris@16:             xa = xa + w(:,k,r) * squeeze(zh(r,k,:))';
Chris@16:         end; 
Chris@16:     end;
Chris@16:     Q = x ./ xa;
Chris@16:     
Chris@16:     % M-step (update h,z)
Chris@16:     if  (lh && lz)
Chris@16:         nh=zeros(K,N);
Chris@16:         for k=1:K
Chris@16:             for r=1:R
Chris@16:                 nh(k,:) = nh(k,:) + squeeze(z(r,k,:))' .* (squeeze(w(:,k,r))' * Q);
Chris@16:                 nz = h(k,:) .* (squeeze(w(:,k,r))' * Q);
Chris@16:                 nz = nz .* squeeze(z(r,k,:))';
Chris@16:                 z(r,k,:) = nz;
Chris@16:             end;
Chris@16:         end;
Chris@16:         nh = h .* nh;
Chris@16:     end
Chris@16:     
Chris@16:     
Chris@16:     % Assign and normalize
Chris@16:     k=1:K;
Chris@16:     n=1:N;
Chris@16:     if lh
Chris@16:         nh = nh.^sh;
Chris@16:         h(:,n) = repmat(sumx(n),K,1) .* (nh(:,n) ./ repmat( sum( nh(:,n), 1), K, 1));
Chris@16:     end
Chris@16:     if lz
Chris@16:         z = z.^sz;
Chris@16:         z(:,k,n) = z(:,k,n) ./ repmat( sum( z(:,k,n), 1), R, 1);
Chris@16:     end
Chris@16:     
Chris@16: end
Chris@16: 
Chris@16: % Show me
Chris@16: if pl
Chris@16:     subplot(3, 1, 1), imagesc(x), axis xy
Chris@16:     subplot(3, 1, 2), imagesc(xa), axis xy
Chris@16:     subplot(3, 1, 3), imagesc(h), axis xy
Chris@16: end