% Example 5-3: Decision algorithm: Which one of two
% non-square input signal matrices (X1 or X2) has been
% transmitted over 2-by-3 MIMO channel?
clear all
K = 4; % Number of time samples (I/O observations)
N = 2; % Number of channel inputs
M = 3; % Number of channel outputs
%====Possible Inputs (a-priori known)===== 
X1 = [4  3
      7  2
      2 -3
      1  5];% (known) input signal matrix number 1
   
X2 = [1  2  
      3 -1  
      7  1  
      2  6];% (known) input signal matrix number 2
%================================   
if size(X1,1) ~= size(X2,1)
   error('Row dimensions of X1 and X2 must be the same.')
end

%=====Output observations========
Y = [1    4   2
     8   -1   1
     9    3   8
     2    7  -3]% Measured output signal matrix  
%================================  
  
[Ux1, ignore] = qr(X1) % QR-factorize input matrix #1
% Ux1 = Orthonormal input basis #1
Wx1 = Ux1*X1; % Generalized Fourier coefficients
%(GFC's) of input #1
Wy1 = Ux1*Y; % Generalized Fourier coefficients
% (GFC's) of output Y in basis Ux1

[Ux2, ignore] = qr(X2) % QR-factorize input matrix #2
% Ux2 = Orthonormal input basis #2
Wx2 = Ux2*X2; % GFC's of input #2
Wy2 = Ux2*Y; % GFC's of output Y in basis Ux2

theta_x = subspace(X1, X2) % angle between subspaces
% spanned by X1 and X2
theta_u = subspace(Ux1,Ux2) % angle between subspaces
% spanned by Ux1 and Ux2

Wc1 = pinv(Wx1)*Wy1; % Channel matrix estimate in basis Ux1
Wc2 = pinv(Wx2)*Wy2; % Channel matrix estimate in basis Ux2

Xh1 = Y*pinv(Wc1); % Estimate of input matrix X w.r.t. basis Ux1
Xh2 = Y*pinv(Wc2); % Estimate of input matrix X w.r.t. basis Ux2

Wx11h = Ux1*Xh1; % GFC's of Xh1 in basis Ux1
Wx12h = Ux2*Xh1; % GFC's of Xh1 in basis Ux2
Wx21h = Ux1*Xh2; % GFC's of Xh2 in basis Ux1
Wx22h = Ux2*Xh2; % GFC's of Xh2 in basis Ux2

Errx11 = norm(X1-Xh1,'fro'); % Calculate Frobenius
% norm of input signal error matrix w.r.t. X1
Errx12 = norm(X2-Xh1,'fro'); % Calculate Frobenius
% norm of input signal error matrix w.r.t. X2
Errx21 = norm(X1-Xh2,'fro'); % Calculate Frobenius
% norm of input signal error matrix w.r.t. X1
Errx22 = norm(X2-Xh2,'fro'); % Calculate Frobenius
% norm of input signal error matrix w.r.t. X2

Errx = [Errx11 Errx12
        Errx21 Errx22]
% Decision matrix based on X estimates
   
% Calculate Frobenius error norms of GFC's   
Err11 = norm(Wx1-Wx11h,'fro'); 
Err12 = norm(Wx2-Wx12h,'fro'); 
Err21 = norm(Wx1-Wx21h,'fro'); 
Err22 = norm(Wx2-Wx22h,'fro');    

Erru = [Err11 Err12
        Err21 Err22]
% Decision matrix based on Wx estimates   
%==========================  
% Find argmin(Erru)
  mini = min(min(Erru)) % Minimum element of decision matrix Erru 
  diago = [Err11-mini Err22-mini];  
  nz = find(diago); % find nonzero elements of vector diago  
  if nz == 1
     p = 2; % Decision in favor of X2
  else
     p = 1; % Decision in favor of X1
  end
  Decision = p % Decision in favor
  %of X1 (p=1) or X2 (p=2), respectively.