function S=sep_l0min(X,A)
% SEP_L0MIN STFT-domain source separation via lp-norm minimization with
% p->0 and with as many active sources as channels, given the mixing matrix
% for instantaneous mixtures or the frequency-dependent mixing matrix for
% convolutive mixtures
%
% S=sep_l0min(X,A)
%
% Inputs:
% X: nbin x nfram x nchan matrix containing mixture STFT coefficients with
% nbin frequency bins, nfram time frames and nchan channels
% A: either a nchan x nsrc mixing matrix (for instantaneous mixtures) or a
% nchan x nsrc x nbin frequency-dependent mixing matrix (for convolutive
% mixtures)
%
% Output:
% S: nbin x nfram x nsrc matrix containing the estimated source STFT
% coefficients
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2008 Emmanuel Vincent
% This software is distributed under the terms of the GNU Public License
% version 3 (http://www.gnu.org/licenses/gpl.txt)
% If you find it useful, please cite the following reference:
% Emmanuel Vincent, "Complex nonconvex lp norm minimization for
% underdetermined source separation," In Proc. Int. Conf. on Independent
% Component Analysis and Blind Source Separation (ICA), pp. 430-437, 2007.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Errors and warnings %%%
if nargin<2, error('Not enough input arguments.'); end
[nbin,nfram,nchan]=size(X);
if nbin==2*floor(nbin/2), error('The number of frequency bins must be odd.'); end
[nchan2,nsrc,nbin2]=size(A);
if nchan2~=nchan, error('The number of channels is different within the mixture and the mixing matrix.'); end
if nbin2==1,
A=repmat(A,[1 1 nbin]);
elseif nbin2~=nbin,
error('The number of frequency bins is different within the mixture and the frequency-dependent mixing matrix.');
end
%%% l0-norm minimization %%%
% Source activity patterns
act=zeros(2^nsrc,nsrc);
for j=1:nsrc,
act(:,j)=reshape(repmat([0 1],[2^(nsrc-j) 2^(j-1)]),2^nsrc,1);
end
act=logical(act(sum(act,2)==nchan,:));
npat=size(act,1);
S=zeros(nbin,nfram,nsrc);
for f=1:nbin,
% Source coefficients
Xf=reshape(X(f,:,:),nfram,nchan).';
SPf=zeros(nchan,nfram,npat);
for p=1:npat,
SPf(:,:,p)=inv(A(:,act(p,:),f))*Xf;
end
% Selection of the best activity pattern
[val,ind]=min(reshape(prod(abs(SPf),1),nfram,npat),[],2);
Sf=zeros(nsrc,nfram);
for p=1:npat,
Sf(act(p,:),ind==p)=SPf(:,ind==p,p);
end
S(f,:,:)=reshape(Sf.',1,nfram,nsrc);
end
return;