acma.m

Analytical constant-modulus algorithm, to separate linear combinations of CM sourcesThe algorithm i

function [A,S,sp] = acma(X,d,d1);
% function [A,S,sp] = acma(X,d,d1);
%
% Analytical constant-modulus algorithm, to separate 
% linear combinations of CM sources
%
% Given data matrix X of rank d, find factorization  
%   W X = S, A = pinv(W), where W,S are rank d1,
% and | S_ij | = 1, for d1 of the signals
%
% d: number of sources;  d1: number of constant-modulus sources
%
% sp: singular values of P (for statistics).  CM sources correspond to zero
% singular values.
%
% See IEEE Tr SP, May 1996

% Alle-Jan van der Veen, Stanford Univ, April 1994
% [e-mail],  [e-mail]


[m,N] = size(X);

if (nargin <= 1), 
   d = m;		% assume fully loaded
end
if (nargin <= 2),
   d1 = d;		% assume all sources are CM
end
flops(0);


% disp('STEP 1: Estimate signal subspace (SVD of X)')
% ---------------------------------------------------------------------------

[u,s,v] = svd(X);		% estimate row space of X
U1 = u(:,1:d);
S1 = s(1:d,1:d);
V1 = v(:,1:d)';			% rows of V1 are basis of S
clear v		% ( v could be large)

% disp(sprintf('flops = %i\n', flops));


% disp('STEP 2: Set up conditions for CM')
% ---------------------------------------------------------------------------

P = zeros(N,d^2);
for i = 1:N,
    v1 = V1(:,i);
    P(i,:) = vecreal(v1*v1')';
end


% Householder transform to map ones(N,1) to e1

v = [1 ; ones(N-1,1)/(1+sqrt(N))];
P = P - 2*v*(v'*P)/(v'*v);	% apply householder transform to P

P = P(2:N,:);			% Throw away first condition
				% We are left with solving Py = 0, for 
				% y = vec(w'*w)



% disp(sprintf('flops = %i\n', flops));


% disp('STEP 3: Solve conditions for CM (SVD of P)')
% ---------------------------------------------------------------------------

% [up,sp,vp] = svd(P);		% solve LS problem P y = 0 (ie find kernel of P)
R = triu(qr(P));		% (qr first, to avoid storing large matrices)
[up,sp,vp] = svd(R);

sp = diag(sp);

% plot(sp,'+');		% CM signals corr. to small singular values


% collect solutions in matrix y
y = zeros(d,d*d1);
for i = 1:d1,
  yi = unvecreal(vp(:,d^2+1-i)); 	% select vectors from the kernel of P
  y(:,(i-1)*d+1:i*d) = yi;	% store in block vector y
end

% At this point, y = [Y_1 ... Y_d1], with Y_k complex, symmetric

% disp(sprintf('flops = %i\n', flops));


% disp('STEP 4: Construct w')
% ---------------------------------------------------------------------------

% Ideally, each yi has rank 1 and is equal to yi = wi'*wi
% However, linear combinations of the yi are also vectors in the kernel
% So, solve lambda_1 y1 + ... + lambda_d yd = wi'*wi (rank 1), for i=1..d.

w = supereig(y);		% decouple solutions to obtain weight vectors

% w: d1 by d

% scale each row of w to norm sqrt(N)
for i=1:d1,
   w(i,:) = w(i,:)/norm(w(i,:))*sqrt(N);
end


% disp(sprintf('flops = %i\n', flops));



% disp('STEP 5: Construct S')
% ---------------------------------------------------------------------------

S = w*V1;

% disp(sprintf('flops = %i\n', flops));

% The remainder is only necc. if we want to know W and A as well
% (note: w was with respect to V, not to X)
W = w*inv(S1)*U1';		% same as Xinv = pinv(X);  W = S*Xinv;
A = pinv(W);			% corresponding array response matrix