PLCA

function [w zf zt h] = coupled_plca(spec_f, spec_t, Bs, Bt, R, iter)
% Perform coupled PLCA
%
% [w zf zt h] = coupled_plca(spec_f, spec_t, Bs, Bt, R, iter)
%
% spec_f: spectrogram with high spectral resolution
% spec_t: spectrogram with high time resolution
% Bs = blurring function for spectral components
% Bt = blurring function for temporal components
% R = number of components
% iter = number of EM iterations
%
% Output parameters
% w = spectral basis vectors (hi-res)
% zf = component weights for spectral vectors
% zt = component weights for temporal vectors
% h = temporal vectors (hi-res too!)
%
% Juhan Nam
%   Apr-11-2010: initial version
%   Mar-11-2013: normalizing blur_h and blur_w right after blurring

% size check
if size(spec_t,1) ~= size(spec_f,1)
    error('coupled_plca.m: the sizes of row are equal');
end

if size(spec_t,2) ~= size(spec_f,2)
    error('coupled_plca.m: the sizes of column are equal');
end

% initializations
[m, n] = size(spec_t);

% spectral basis vectors,
%w = rand(m, R)+1;
ind = randperm(n);
w = spec_f(:,ind(1:R))+eps;
w = bsxfun(@rdivide, w, sum( w, 1));

% temporal basis vectors
h = rand(R, n)+1;
h = bsxfun(@rdivide, h, sum( h, 2));

zf = rand(R)+1; % weights vector for good freq res plca
zf = zf / sum(zf); % normalize

zt = rand(R)+1; % weights vector for good time res plca
zt = zt / sum(zt); % normalize


Bs = Bs(:);  % vertical
Bt = Bt(:)'; % horizontal

% EM iterations
for it = 1:iter
    % E-step for high frequency resolution spectrogram
    blur_h = conv2(h, Bt, 'same');
    blur_h = bsxfun(@rdivide, blur_h, sum(blur_h,2)+eps);
    zfh = diag(zf) * blur_h;
    Rf = spec_f ./ (w * zfh);

    % E for high time resolution spectrogram
    zth = diag(zt) * h;
    blur_w = conv2(w, Bs, 'same');
    blur_w = bsxfun(@rdivide, blur_w, sum(blur_w,1)+eps);
    Rt = spec_t ./ (blur_w * zth);

    % M-step
    nw = w .* ( Rf * zfh');
    nh = zth .* (w'* Rt);
    nzf = sum(nw, 1);
    nzt = sum( blur_w, 1);

    % normalize
    w = bsxfun(@rdivide, nw, sum(nw,1)+eps);
    h = bsxfun(@rdivide, nh, sum(nh,2)+eps);
    zf = nzf / sum(nzf+eps);
    zt = nzt / sum(nzt+eps);
end