% ========================================================================
%
% A simple offline implementation of
%
% Unconstrained Multivariate Frequency-Domain ICA for Convolutive Mixtures
%
% Herbert Buchner, Apr. 2011                           www.buchner-net.com
%
% ========================================================================

%-------------------------------------------------------
%-------------parameters for the algorithm--------------
%-------------------------------------------------------
cfg.nmic = 2;        % number of microphones
cfg.nsrc = 2;        % number of sources
%cfg.fs = 16000;      % sampling frequency of the microphone signals

cfg.fftsize = [1024 512];  % frameSize, shift
cfg.bins = cfg.fftsize(1)/2+1;

cfg.loops = 250;     % number of iterations for ICA part
cfg.eta = 0.01;      % stepsize
cfg.delta = eps;     % regularization for multivariate score function

%---------------------------------------------------------------------
%----------------load microphone signals x----------------------------
%---------------------------------------------------------------------
fprintf('reading inputs...\n')
for p=1:cfg.nmic
  [x(:,p), cfg.fs, tmp] = wavread(strcat('x', num2str(p), '.wav'));
end

%---------------------------------------------------------------------
%------------calculate frequency-domain version of x------------------
%---------------------------------------------------------------------
fprintf('STFT:   frameSize = %d, shift = %d\n', cfg.fftsize(1), cfg.fftsize(2));
hwindow = hanning(cfg.fftsize(1), 'periodic');

for p=1:cfg.nmic
   tmp = sfft(x(:,p), hwindow, cfg.fftsize(2));
   X(p,:,:) = shiftdim( tmp(1:cfg.bins,:), 1 );   % X = zeros(cfg.nmic, nFrame, cfg.bins);
end
nFrame = size(X,2);  % number of frames

% normalization of microphone signals
for p = 1:cfg.nmic
   for f =1:cfg.bins
       X(p,:,f) = X(p,:,f)- mean(squeeze(X(p,:,f)));
   end
end

%---------------------------------------------------------------------
%--------------------Compute demixing filter W------------------------
%---------------------------------------------------------------------
W = zeros(cfg.bins, cfg.nsrc, cfg.nmic);
Y = zeros(cfg.nsrc, nFrame, cfg.bins);
Yf = [];

%---initialization---
Wt = zeros(cfg.fftsize(1),cfg.nsrc,cfg.nmic);
for p = 1:cfg.nmic
   Wt(1,p,p) = 1;
end
W = fft(Wt,[],1);


%------------main loop----------------
for i=1:cfg.loops
   fprintf('Estimation:    iteration %d of %d\nM', i, cfg.loops);

   for f=1:cfg.bins,
       Wf = squeeze(W(f,:,:));
       Xf = X(:,:,f);
       Yf = Wf * Xf;    % circular convolution to obtain output signals
       Y(:,:,f) = Yf;
   end

   % !!!!!!! norm necessary for MULTIvariate pdf !!!!!!! 
   phiYfvec = sqrt(sum((abs(Y).^2),3));

   for f=1:cfg.bins
       Wf = squeeze(W(f,:,:));
       Xf = X(:,:,f);
       Yf = Y(:,:,f);

       % !!!!!!! nonlinearity based on UNIvariate pdf !!!!!!! 
       phiYf = sqrt(cfg.bins)*Yf./(abs(Yf)+cfg.delta);

       % !!!!!!! nonlinearity based on MULTIvariate pdf !!!!!!! 
%        phiYf = sqrt(cfg.bins)*Yf./(phiYfvec+cfg.delta);

       %-------------------------------------
       % natural gradient (holonomic version)
       %-------------------------------------
       Syy_norm = (phiYf * Yf') ./ nFrame;
       Wf = Wf + cfg.eta * ( eye(cfg.nsrc)-Syy_norm ) * Wf;

       %-------------------------------------------------
       % Minimal distortion principle (Matsuoka, ICA2001)
       %-------------------------------------------------
       Wfinv = pinv(Wf);
       W(f,:,:) = diag(diag(Wfinv)) * Wf;
   end
end
fprintf('\n');

%---------------------------------------------------------------------
%----------Compute FIR demixing filters from freq-domain W---------------------
%---------------------------------------------------------------------
Wt = real(ifft( [W(1:cfg.bins,:,:); conj(W(cfg.bins-1:-1:2,:,:))], cfg.fftsize(1)));
Wt = circshift(Wt,[cfg.fftsize(1)/2,0,0]);  % put peaks in the middle of filter

%---------------------------------------------------------------------
%---------Compute output signals--------------------------------------
%---------------------------------------------------------------------
for p=1:cfg.nsrc,
 % make initial y(:,p) with the first microphone
 y(:,p) = fftfilt(Wt(:,p,1), x(:,1));
 % loop for the second or more.  Assume the lengths of x{q} are the same.
 for q=2:cfg.nmic,
   y(:,p) = y(:,p) + fftfilt(Wt(:,p,q), x(:,q));
 end
 % save output signals
 wavwrite(y(:,p)./(1.01*max(abs(y(:,p)))), cfg.fs, strcat('y', num2str(p), '.wav'))
end