function [y,w] = ECoBliSS(xn,N,alpha,B,P,L2)
% ECoBliSS - Adaptive Blind Source Separation
%
% usage: [y,w,e] = ECoBliSS(xn,N,alpha,B,P,L2)
%
% xn - sensor data
% N - #tabs per sensor (should be a power of 2)
%alpha - forgetting factor for autocorrelation matrix update (0.9<alpha<1)
% B - block length - Ignored for now, B=N;
% P - weights are updated at every P-th block of data
% L2 - number of blind channels (the first L2 columns of xn represent mixed signals, the remaining columns contain signals the are independent of the signals to be recovered.)
%
% y - array output
% w - final weights
% e - instantaneous error
%
% See also: ECoBliSS_update, RF2RFsq, ECoBliSS_TDwindow, MIRROR, sqrtmH
%
% Author: D.W.E. Schobben (ds@altavista.net)
% initialisation
y = zeros(length(xn),L2);
J = min(size(xn)); % #inputs
M = 2*N; %%% 2^ceil(log2(N+B-1)); % transform length M=nearest power of 2
W = zeros(M,L2,J);
RF = zeros(M,J,J);
for ii = 1:J; RF(:,ii,ii) = ones(M,1); end
imax = length(y)-M+1; % #iterations
for n=1:J, xn(:,n) = xn(:,n) - mean(xn(:,n)); end
% FD-MC-BSS
for i = 1:B:imax,
% Transform data to frequency domain
xx = xn(i:i+M-1,:); % sensor data starting with oldest samples
XX = fft(xx); % frequency domain microphone data
% Filter data in FD
if i>5,
for ii=1:L2,
temp = real(ifft( sum((XX.*permute(W(:,ii,:),[1 3 2]))')' ));
y(i+M-B:i+M-1,ii) = temp(M-B+1:M);
end
end
% Update covariance matrix in the frequency domain
for r = 1:J, for l = 1:J,
RF(:,r,l) = alpha*RF(:,r,l) + (1-alpha)*( conj(XX(:,r)).*XX(:,l).*(B/M) );
end; end
if ( (rem((i-1)/B,P) == 1 | P==1) );
fprintf('Percentage completed: %0.2g\n',100*i/imax) % Print progress to stdout
% Decompose RF
RFsq = RF2RFsq(permute(RF,[2 3 1])); % Pretty expensive (1.28s)
% Update the FD weights
W = ECoBliSS_update(W,RFsq,N,RF); % Most expensive (2.09s)
end
end
fprintf('Percentage completed: done\n');
% Calculate TD weights
for r=1:L2, for s=1:J, temp = real(sqrt(M)*ifft(W(:,r,s))); w(:,r,s) = flipud(temp(1:N)); end; end
