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% allejan@isl.stanford.edu,  allejan@cas.et.tudelft.nl[m,N] = size(X);if (nargin <= 1),    d = m;		% assume fully loadedendif (nargin <= 2),   d1 = d;		% assume all sources are CMendflops(0);% disp('STEP 1: Estimate signal subspace (SVD of X)')% ---------------------------------------------------------------------------[u,s,v] = svd(X);		% estimate row space of XU1 = u(:,1:d);S1 = s(1:d,1:d);V1 = v(:,1:d)';			% rows of V1 are basis of Sclear 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 e1v = [1 ; ones(N-1,1)/(1+sqrt(N))];P = P - 2*v*(v'*P)/(v'*v);	% apply householder transform to PP = 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 yy = 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 yend% 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