www.pudn.com > convolutive.rar > convbss.m

function [y,Wt,J] = convbss(x,T,Q,K,N,verbose,iter,printevery,eta,ds); 
[rx,cx] = size(x); 
 
% assign default values to the unassigned 
if ~exist('T')|isempty(T), T = 512; end 
if mod(T,2), error('T must be multiple of 2'), end 
hT = T/2+1; % the relevant part of the FFT 
if ~exist('Q')|isempty(Q), Q = 128; end 
if ~exist('K')|isempty(K), K = 5; end 
if ~exist('N')|isempty(N), N = floor(rx/T/K); end 
if ~exist('verbose')|isempty(verbose), verbose = 1; end 
if ~exist('iter')|isempty(iter), iter = 1000; end 
if ~exist('printevery')|isempty(printevery), printevery = 10; end 
if ~exist('eta')|isempty(eta), eta = 1.0; end 
if ~exist('ds')|isempty(ds), ds = cx; end 
 
% we want some of the large variables to be global 
global Rx Rxk m dsdshT dsdxhT 
 
% the mixed signals must have zero mean and variance one 
x=(x-repmat(mean(x),length(x),1))./std(x(:));  
 
%x = (x-repmat(mean(x),[rx 1]))./repmat(std(x),[rx 1]); 
dx = cx; % the number of senors is the number of columns 
if ds>dx, error('ds is not allowed to be greater than dx'), end 
 
% calculate the cross-power-spectrum Rx 
Rx = cps(x,T,hT,K,N,dx,verbose); % size(Rx) == [dx dx hT K] 
 
% calculate some straightforward power normalization 
Rxk = sum(Rx,4); % sum over the k matrices 
m = 1./repmat(sum(sum(Rxk.*conj(Rxk),2),1),[ds dx 1]); 
 
% initialize the demixing filter to the identity filter 
Wt = zeros(ds,dx,T); % Wt = W in time domain 
Wt(:,:,1) = eye(ds,dx); 
Wf = shortfft(Wt); % Wf = W in frequency domain 
 
% define the indices of the diagonal of some matrices 
dsdshT = find(repmat(eye(ds,ds),[1 1 hT])); 
dsdxhT = find(repmat(eye(ds,dx),[1 1 hT])); 
 
oldJ = inf; 
 
% do the gradient descent 
while iter > 0 
if verbose, fprintf('iter=%d \r',iter), end 
 
if verbose & (mod(iter,printevery)==0 | iter==1) 
% calculate the error 
J = calcJ(Wf,hT,K); 
 
% check whether we still made progress 
if oldJ % we are no longer making progress, let's break this loop  
break 
else 
oldJ = J; 
end 
fprintf('iter=%d J=%f \n',iter,J) 
end 
 
% calculate the gradient of the error 
dWf = calcdWf(Wf,ds,dx,hT,K); 
 
% update the unmixing filter Wf 
Wf = Wf - eta*dWf; 
 
% project the filter matrix such that Wt(:,:,Q:end) == 0 
Wt = shortifft(Wf); 
Wt(:,:,Q+1:end) = 0; 
Wf = shortfft(Wt); 
 
iter = iter - 1; 
end 
 
% calculate the error 
J = calcJ(Wf,hT,K); 
 
% finally calculate the demixed signals 
if verbose, fprintf('calculate the demixed signals (this might take a while)...'), end 
Wt = shortifft(Wf); 
y = demix(Wt,x); 
if verbose, fprintf('done\n'), end