1 % ========================================================================
   2 %
   3 % A simple offline implementation of
   4 %
   5 % Unconstrained Multivariate Frequency-Domain ICA for Convolutive Mixtures
   6 %
   7 % Herbert Buchner, Apr. 2011                           www.buchner-net.com
   8 %
   9 % ========================================================================
  10 
  11 %-------------------------------------------------------
  12 %-------------parameters for the algorithm--------------
  13 %-------------------------------------------------------
  14 cfg.nmic = 2;        % number of microphones
  15 cfg.nsrc = 2;        % number of sources
  16 %cfg.fs = 16000;      % sampling frequency of the microphone signals
  17 
  18 cfg.fftsize = [1024 512];  % frameSize, shift
  19 cfg.bins = cfg.fftsize(1)/2+1;
  20 
  21 cfg.loops = 250;     % number of iterations for ICA part
  22 cfg.eta = 0.01;      % stepsize
  23 cfg.delta = eps;     % regularization for multivariate score function
  24 
  25 %---------------------------------------------------------------------
  26 %----------------load microphone signals x----------------------------
  27 %---------------------------------------------------------------------
  28 fprintf('reading inputs...\n')
  29 for p=1:cfg.nmic
  30   [x(:,p), cfg.fs, tmp] = wavread(strcat('x', num2str(p), '.wav'));
  31 end
  32 
  33 %---------------------------------------------------------------------
  34 %------------calculate frequency-domain version of x------------------
  35 %---------------------------------------------------------------------
  36 fprintf('STFT:   frameSize = %d, shift = %d\n', cfg.fftsize(1), cfg.fftsize(2));
  37 hwindow = hanning(cfg.fftsize(1), 'periodic');
  38 
  39 for p=1:cfg.nmic
  40    tmp = sfft(x(:,p), hwindow, cfg.fftsize(2));
  41    X(p,:,:) = shiftdim( tmp(1:cfg.bins,:), 1 );   % X = zeros(cfg.nmic, nFrame, cfg.bins);
  42 end
  43 nFrame = size(X,2);  % number of frames
  44 
  45 % normalization of microphone signals
  46 for p = 1:cfg.nmic
  47    for f =1:cfg.bins
  48        X(p,:,f) = X(p,:,f)- mean(squeeze(X(p,:,f)));
  49    end
  50 end
  51 
  52 %---------------------------------------------------------------------
  53 %--------------------Compute demixing filter W------------------------
  54 %---------------------------------------------------------------------
  55 W = zeros(cfg.bins, cfg.nsrc, cfg.nmic);
  56 Y = zeros(cfg.nsrc, nFrame, cfg.bins);
  57 Yf = [];
  58 
  59 %---initialization---
  60 Wt = zeros(cfg.fftsize(1),cfg.nsrc,cfg.nmic);
  61 for p = 1:cfg.nmic
  62    Wt(1,p,p) = 1;
  63 end
  64 W = fft(Wt,[],1);
  65 
  66 
  67 %------------main loop----------------
  68 for i=1:cfg.loops
  69    fprintf('Estimation:    iteration %d of %d\nM', i, cfg.loops);
  70 
  71    for f=1:cfg.bins,
  72        Wf = squeeze(W(f,:,:));
  73        Xf = X(:,:,f);
  74        Yf = Wf * Xf;    % circular convolution to obtain output signals
  75        Y(:,:,f) = Yf;
  76    end
  77 
  78    % !!!!!!! norm necessary for MULTIvariate pdf !!!!!!! 
  79    phiYfvec = sqrt(sum((abs(Y).^2),3));
  80 
  81    for f=1:cfg.bins
  82        Wf = squeeze(W(f,:,:));
  83        Xf = X(:,:,f);
  84        Yf = Y(:,:,f);
  85 
  86        % !!!!!!! nonlinearity based on UNIvariate pdf !!!!!!! 
  87        phiYf = sqrt(cfg.bins)*Yf./(abs(Yf)+cfg.delta);
  88 
  89        % !!!!!!! nonlinearity based on MULTIvariate pdf !!!!!!! 
  90 %        phiYf = sqrt(cfg.bins)*Yf./(phiYfvec+cfg.delta);
  91 
  92        %-------------------------------------
  93        % natural gradient (holonomic version)
  94        %-------------------------------------
  95        Syy_norm = (phiYf * Yf') ./ nFrame;
  96        Wf = Wf + cfg.eta * ( eye(cfg.nsrc)-Syy_norm ) * Wf;
  97 
  98        %-------------------------------------------------
  99        % Minimal distortion principle (Matsuoka, ICA2001)
 100        %-------------------------------------------------
 101        Wfinv = pinv(Wf);
 102        W(f,:,:) = diag(diag(Wfinv)) * Wf;
 103    end
 104 end
 105 fprintf('\n');
 106 
 107 %---------------------------------------------------------------------
 108 %----------Compute FIR demixing filters from freq-domain W---------------------
 109 %---------------------------------------------------------------------
 110 Wt = real(ifft( [W(1:cfg.bins,:,:); conj(W(cfg.bins-1:-1:2,:,:))], cfg.fftsize(1)));
 111 Wt = circshift(Wt,[cfg.fftsize(1)/2,0,0]);  % put peaks in the middle of filter
 112 
 113 %---------------------------------------------------------------------
 114 %---------Compute output signals--------------------------------------
 115 %---------------------------------------------------------------------
 116 for p=1:cfg.nsrc,
 117  % make initial y(:,p) with the first microphone
 118  y(:,p) = fftfilt(Wt(:,p,1), x(:,1));
 119  % loop for the second or more.  Assume the lengths of x{q} are the same.
 120  for q=2:cfg.nmic,
 121    y(:,p) = y(:,p) + fftfilt(Wt(:,p,q), x(:,q));
 122  end
 123  % save output signals
 124  wavwrite(y(:,p)./(1.01*max(abs(y(:,p)))), cfg.fs, strcat('y', num2str(p), '.wav'))
 125 end