convbss.asv

本程序自己谢的卷积盲分离的程序

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]
%加了
%function Rx = cps(x,T,hT,K,N,dx,verbose);
% calculate the cross-power-spectrum Rx

% check whether these parameters can be used with signals of length rx,
% this check is done by checking whether the last index that will be
% accessed during the calculation of Rx exists
if K*T*N > size(x,1)
error('the signals are too short for the chosen parameters')
end

% do the calculation
if verbose, disp('estimate the cross-power spectrum of x'), end
Rx = zeros(dx,dx,hT,K);
XX = zeros(dx,dx,hT);
for k=1:K
for n=1:N
if verbose, fprintf('%d/%d %d/%d \r',k,K,n,N); end
indx = (k-1)*T*N + (n-1)*T + (1:T);
X=fft(x(indx,:));
% since we are only interested in filter in frequency domain
% that correspond to real valued filters in time domain,
% we can ignore the negative frequencies
X=X(1:hT,:);
for row = 1:dx, for col = 1:dx
XX(row,col,:) = reshape(conj(X(:,col)).*X(:,row),[1 1 hT]);
end, end
Rx(:,:,:,k) = Rx(:,:,:,k) + XX;
end;
end;
Rx = Rx/N/T;



% 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
%加了
%function Wt = shortifft(Wf)
Wf = real(ifft(cat(3,Wt,conj(Wt(:,:,[end-1:-1:2]))),[],3));



% 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);
%function J = calcJ(Wf,hT,K);
global Rx dsdshT
% calculate the error term for diagonalization
J = 0;
for k = 1:K
%加了
%function C = fast_ip(A,B)
% for high-(more than three)-dimensional arrays with a lot
% of elements in those higher dimensions it might be very 
% expensive to calculate the inner product. this function 
% implements it in a fast way
[rA cA sA] = size(Wf);
[rB cB sB] = size(Rx(:,:,:,k);
if cA~=rB, error('the inner dimensions of A and B do not agree'), end
if sA~=sB, error('the higher dimensions of A and B do not agree'), end
C = zeros(rA,cB,sB);
for row=1:rA, for col=1:cB
C(row,col,:) = sum(reshape(A(row,:,:),[cA 1 sA]) .* B(:,col,:),1);
end, end
%Eomega = fast_ip(fast_ip(Wf,Rx(:,:,:,k)),hermite3(Wf));

%function Ah = hermite3(A)
% for three dimensional arrays the function ' can not be used
% to conjugate-transpose (hermite) a matrix
Q = permute(conj(Wf),[2 1 3]);
%function C = fast_ip(A,B)
% for high-(more than three)-dimensional arrays with a lot
% of elements in those higher dimensions it might be very 
% expensive to calculate the inner product. this function 
% implements it in a fast way
[rA cA sA] = size(C);
[rB cB sB] = size(Q);
if cA~=rB, error('the inner dimensions of A and B do not agree'), end
if sA~=sB, error('the higher dimensions of A and B do not agree'), end
Eomega = zeros(rA,cB,sB);
for row=1:rA, for col=1:cB
Eomega(row,col,:) = sum(reshape(A(row,:,:),[cA 1 sA]) .* B(:,col,:),1);
end, end
Eomega(dsdshT) = 0; % this uses the optimal lambdaS
% note that trace(A*A') == A(:)' * A(:)
J = J + Eomega(:)' * Eomega(:);
end
J = real(J)/K/hT;


% 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);
%function dWf = calcdWf(Wf,ds,dx,hT,K);
%global Rx m dsdshT dsdxhT

% calculate the gradient of the objective function
dWf = zeros(ds,dx,hT);
%加了
%hWf = hermite3(Wf);
%function Ah = hermite3(A)
% for three dimensional arrays the function ' can not be used
% to conjugate-transpose (hermite) a matrix
hWf = permute(conj(Wf),[2 1 3]);


for k = 1:K
WfRx = fast_ip(Wf,Rx(:,:,:,k));
Eomega = fast_ip(WfRx,hWf);
Eomega(dsdshT) = 0; % this uses the optimal lambdaS
dWf = dWf + fast_ip(Eomega,WfRx);
end
dWf = 2*dWf;

%% apply the straightforward power normalization
dWf = m.*dWf;

% enforce unit gains of the diagonal filters
dWf(dsdxhT) = 0;



% update the unmixing filter Wf
Wf = Wf - eta*dWf;

% project the filter matrix such that Wt(:,:,Q:end) == 0
%加了
%Wt = shortifft(Wf);
%function Wt = shortifft(Wf)
Wt = real(ifft(cat(3,Wf,conj(Wf(:,:,[end-1:-1:2]))),[],3));


Wt(:,:,Q+1:end) = 0;
%加了
%Wf = shortfft(Wt);
%function Wf = shortfft(Wt)
Wf = fft(Wt,[],3);
Wf = Wf(:,:,1:end/2+1);

iter = iter - 1;
end

% calculate the error
%加了
%J = calcJ(Wf,hT,K);
%function J = calcJ(Wf,hT,K);
%global Rx dsdshT
% calculate the error term for diagonalization
J = 0;
for k = 1:K
Eomega = fast_ip(fast_ip(Wf,Rx(:,:,:,k)),hermite3(Wf));
Eomega(dsdshT) = 0; % this uses the optimal lambdaS
% note that trace(A*A') == A(:)' * A(:)
J = J + Eomega(:)' * Eomega(:);
end
J = real(J)/K/hT;

% finally calculate the demixed signals
if verbose, fprintf('calculate the demixed signals (this might take a while)...'), end
%加了
%Wt = shortifft(Wf);
%function Wt = shortifft(Wf)
Wt = real(ifft(cat(3,Wf,conj(Wf(:,:,[end-1:-1:2]))),[],3));


%加了
%y = demix(Wt,x);
%function y = demix(Wt,x)
% this is very slow, better would be something 
% like overlap-save, see Haykin or Oppenheim/Schaefer
% it is assumed that size(Wt) = [ds dx T]
% size(x) = [?? ds]

[ds dx T] = size(Wt);
[rx cx] = size(x);

% check for consistency
if ds ~= cx
error('Wt and x must fit together')
end

% allocate memory for y
y = zeros(rx,ds);

% do the convolution using matlab's built-in filter
for ii = 1:ds
for jj = 1:dx
y(:,ii) = y(:,ii) + filter(Wt(ii,jj,:),1,x(:,jj));
end
end

if verbose, fprintf('done\n'), end