Chris@2: function [w,h,z,u,xa] = cplcaMT( x, K, T, R, w, h, z, u, iter, sw, sh, sz, su, lw, lh, lz, lu, pa) Chris@2: % function [w,h,xa2] = cplcaMT( x, K, T, R, w, h, z, u, iter, sw, sh, sz, su, lw, lh, lz, lu) Chris@2: % Chris@2: % Perform multiple-source, multiple-template SIPLCA for transcription Chris@2: % Chris@2: % Inputs: Chris@2: % x input distribution Chris@2: % K number of components Chris@2: % T size of components Chris@2: % R size of sources Chris@2: % w initial value of p(w) [default = random] Chris@2: % h initial value of p(h) [default = random] Chris@2: % z initial value of p(z) [default = random] Chris@2: % iter number of EM iterations [default = 10] Chris@2: % sw sparsity parameter for w [default = 1] Chris@2: % sh sparsity parameter for h [default = 1] Chris@2: % sz sparsity parameter for z [default = 1] Chris@2: % lw flag to update w [default = 1] Chris@2: % lh flag to update h [default = 1] Chris@2: % lh flag to update h [default = 1] Chris@2: % pa source-component activity range [Rx2] Chris@2: % Chris@2: % Outputs: Chris@2: % w p(w) - spectral bases Chris@2: % h p(h) - pitch impulse Chris@2: % z p(z) - mixing matrix for p(h) Chris@2: % xa approximation of input Chris@2: Chris@2: % Emmanouil Benetos 2011, based on cplca code by Paris Smaragdis Chris@2: Chris@2: Chris@2: % Sort out the sizes Chris@2: wc = 2*size(x)-T; Chris@2: hc = size(x)+T-1; Chris@2: Chris@2: % Default training iterations Chris@2: if ~exist( 'iter') Chris@2: iter = 10; Chris@2: end Chris@2: Chris@2: Chris@2: % Initialize Chris@2: sumx = sum(x); Chris@2: Chris@2: if ~exist( 'w') || isempty( w) Chris@2: w = cell(R, K); Chris@2: for k = 1:K Chris@2: for r=1:R Chris@2: w{r,k} = rand( T); Chris@2: w{r,k} = w{r,k} / sum( w{r,k}(:)); Chris@2: end Chris@2: end Chris@2: end Chris@2: if ~exist( 'h') || isempty( h) Chris@2: h = cell(1, K); Chris@2: for k = 1:K Chris@2: h{k} = rand( size(x)-T+1); Chris@2: h{k} = h{k}; Chris@2: end Chris@2: end Chris@2: if ~exist( 'z') || isempty( z) Chris@2: z = cell(1, K); Chris@2: for k = 1:K Chris@2: z{k} = rand( size(x,2),1); Chris@2: z{k} = z{k}; Chris@2: end Chris@2: end Chris@2: if ~exist( 'u') || isempty( u) Chris@2: u = cell(R, K); Chris@2: for k = 1:K Chris@2: for r=1:R Chris@2: if( (pa(r,1) <= k && k <= pa(r,2)) ) Chris@2: u{r,k} = ones( size(x,2),1); Chris@2: else Chris@2: u{r,k} = zeros( size(x,2),1); Chris@2: end Chris@2: end; Chris@2: end Chris@2: end Chris@2: Chris@2: fh = cell(1, K); Chris@2: fw = cell(R, K); Chris@2: for k = 1:K Chris@2: fh{k} = ones(wc) + 1i*ones(wc); Chris@2: for r=1:R Chris@2: fw{r,k} = ones(wc) + 1i*ones(wc); Chris@2: end; Chris@2: end; Chris@2: Chris@2: Chris@2: Chris@2: % Make commands for subsequent multidim operations and initialize fw Chris@2: fnh = 'c(hc(1)-(T(1)+(1:size(h{k},1))-2),hc(2)-(T(2)+(1:size(h{k},2))-2))'; Chris@2: xai = 'xa(1:size(x,1),1:size(x,2))'; Chris@2: flz = 'xbar(end:-1:1,end:-1:1)'; Chris@2: Chris@2: for k = 1:K Chris@2: for r=1:R Chris@2: if( (pa(r,1) <= k && k <= pa(r,2)) ) Chris@2: fw{r,k} = fftn( w{r,k}, wc); Chris@2: end; Chris@2: end; Chris@2: end; Chris@2: Chris@2: % Iterate Chris@2: for it = 1:iter Chris@2: Chris@2: %disp(['Iteration: ' num2str(it)]); Chris@2: Chris@2: % E-step Chris@2: xa = eps; Chris@2: for k = 16:73 Chris@2: fh{k} = fftn( h{k}, wc); Chris@2: for r=1:R Chris@2: if( (pa(r,1) <= k && k <= pa(r,2)) ) Chris@2: xa1 = abs( real( ifftn( fw{r,k} .* fh{k}))); Chris@2: xa = xa + xa1(1:size(x,1),1:size(x,2)) .*repmat(z{k},1,size(x,1))'.*repmat(u{r,k},1,size(x,1))'; Chris@2: clear xa1; Chris@2: end Chris@2: end Chris@2: end Chris@2: Chris@2: xbar = x ./ xa; Chris@2: xbar = eval( flz); Chris@2: fx = fftn( xbar, wc); Chris@2: Chris@2: Chris@2: % M-step Chris@2: for k = 16:73 Chris@2: Chris@2: Chris@2: % Update h, z, u Chris@2: nh=eps; Chris@2: for r=1:R Chris@2: if( (pa(r,1) <= k && k <= pa(r,2)) ) Chris@2: c = abs( real( ifftn( fx .* fw{r,k} ))); Chris@2: nh1 = eval( fnh); Chris@2: nh1 = nh1 .*repmat(u{r,k},1,size(h{k},1))'; Chris@2: nh = nh + nh1; Chris@2: Chris@2: nhu = eval( fnh); Chris@2: nhu = nhu .* h{k}; Chris@2: nu = sum(nhu)'; Chris@2: nu = u{r,k} .* nu + eps; Chris@2: if lu Chris@2: u{r,k} = nu; Chris@2: end; Chris@2: Chris@2: end; Chris@2: end Chris@2: nh = h{k} .* (nh.^sh); Chris@2: nz = sum(nh)'; Chris@2: nz = z{k} .* nz + eps; Chris@2: Chris@2: Chris@2: % Assign and normalize Chris@2: if lh Chris@2: h{k} = nh; Chris@2: end Chris@2: if lz Chris@2: z{k} = nz; Chris@2: end Chris@2: Chris@2: Chris@2: end Chris@2: Chris@2: % Normalize z over t Chris@2: if lz Chris@2: Z=[]; for i=1:K Z=[Z z{i}]; end; Chris@2: Z = Z.^sz; Chris@2: Z(1:end,1:15)=0; Chris@2: Z(1:end,74:88)=0; Chris@2: Z = Z./repmat(sum(Z,2),1,K); z = num2cell(Z,1); %figure; imagesc(imrotate(Z,90)); Chris@2: end Chris@2: Chris@2: % Normalize u over z,t Chris@2: if lu Chris@2: U=[]; for r=1:R U(r,:,:) = cell2mat(u(r,:)); end; Chris@2: for i=1:size(U,2) for j=1:size(U,3) U(:,i,j) = U(:,i,j).^su; U(:,i,j) = U(:,i,j) ./ sum(U(:,i,j)); end; end; Chris@2: U0 = permute(U,[2 1 3]); u = squeeze(num2cell(U0,1)); Chris@2: end Chris@2: Chris@2: % Normalize h over z,t Chris@2: H=[]; for k=1:K H(k,:,:) = cell2mat(h(k)); end; H0 = permute(H,[2 1 3]); Chris@2: for i=1:size(H0,2) for j=1:size(H0,3) H0(:,i,j) = sumx(j)* (H0(:,i,j) ./ sum(H0(:,i,j))); end; end; Chris@2: h = squeeze(num2cell(squeeze(H0),[1 3])); for k=1:K h{k} = squeeze(h{k}); end; Chris@2: Chris@2: %figure; imagesc(imrotate(xa',90)); Chris@2: Chris@2: end Chris@2: Chris@2: %figure; imagesc(imrotate(xa',90));