% A simple CDMA simulation to demo comfac.m

% Reference:
% N. D. Sidiropoulos, G. B. Giannakis, and Rasmus Bro, 
% ``Blind PARAFAC Receivers for DS-CDMA Systems'', 
% IEEE Transactions on Signal Processing, 48(3):810-823, Mar. 2000

% demo by N. Sidiropoulos and R. Bro

clear all;

F=5;  % users
I=4;  % antennas
J=50; % symbol snapshots
K=4;  % spreading gain

A = randn(I,F)+j*randn(I,F); % array response
B = sign(randn(J,F)); % BPSK symbol matrix
for f=1:F, % for scaling ambiguity removal later:
 B(:,f) = B(:,f) / B(1,f);
end
C = randn(K,F); % randomly drawn symbol-periodic signatures

X = zeros(I,J,K);
for k=1:K,
 X(:,:,k) = A*diag(C(k,:))*B.';
end

% add noise if you wish:
% sig2 = 0.01;
% W = (sig2/2)*randn(size(X))+j*(sig2/2)*randn(size(X));
% X = X + W;

%[estA,estC,estB]=TALS(X,F); % plain-vanilla
[estA,estB,estC]=comfac(X,F); % faster and better

% remove permutation and scale ambiguity:

for f=1:F, % scaling ambiguity removal
 estB(:,f) = estB(:,f) / estB(1,f);
end

estB = sign(real(estB)); % hard decision detector

% now remove permutation ambiguity by greedy-LS column matching:

% Calculate permutation that best matches the signal streams
% This is done for performance evaluation purposes only.
% Here we do it in a greedy (suboptimal) way, which is, nevertheless,
% robust & computationally efficient for any number of users.
% define a distance matrix DISTMAT whose element DISTMAT(i,j) is the
% distance between i-th transmitted signal and j-th received signal

DISTMAT=ones(F,F);
for sd=1:F
 for rs=1:F
  DISTMAT(sd,rs)=norm(B(:,sd)-estB(:,rs),'fro');
 end
end
maxD=max(max(DISTMAT))+1;

fS=ones(size(B));
for k=1:F
 row_ps=0;col_ps=0;
 minD=min(min(DISTMAT));
 % the following can be more efficiently done, but what the heck:
 for row=1:F
  for col=1:F
   if (minD==DISTMAT(row,col)) row_ps=row; col_ps=col; end
  end
 end
 fS(:,row_ps)=estB(:,col_ps);
 DISTMAT(row_ps,:)=maxD*ones(size(DISTMAT(row_ps,:)));
 DISTMAT(:,col_ps)=maxD*ones(size(DISTMAT(:,col_ps)));
end

estB = fS;

B
estB

NumBitErrors  = norm(B ~= estB,'fro')^2

BER = NumBitErrors/(J*F)

figure(1);
clf;
for f=1:F,
 subplot(F,1,f);
 bar(B(:,f));
end

figure(2);
clf;
for f=1:F,
 subplot(F,1,f);
 bar(estB(:,f),'r');
end

disp('Blue - true user signals; Red - blind PARAFAC estimates');