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