%% This code snipet is for self interference comparison between the Softnull and MHBFD

H(:,:)=internal_ip2(10,33:64,1:32,9);
[U S V] = svd(H);
h_down(:,:) = downlink(10,:,1:32,9);
%prompt = "What is the number of effective antenna? "
x = 16; %input(prompt);
ef = 32 - x + 1;
p_self = V(:,ef:32); % Softnull precoding
Heff = h_down*p_self;
pdown = Heff'*((Heff*Heff')^(-1));  % closed form formula for zero-forcing
s_down = ones(4,1);  % sending four "1" symbol for each client
MP =(p_self*pdown*s_down);
x_down = (1/(norm(MP,'fro')^2))*MP;
SI_term = H*x_down; 
SI_power = db(norm(SI_term,'fro')^2);


%% 
i=10;
h1(:,:) = internal_ip2(i,33:64,1:8,9);
h2(:,:) = internal_ip2(i,33:64,9:16,9);
h3(:,:) = internal_ip2(i,33:64,17:24,9);
h4(:,:) = internal_ip2(i,33:64,25:32,9);
[U1 S1 V1] = svd(h1);
[U2 S2 V2] = svd(h2);
[U3 S3 V3] = svd(h3);
[U4 S4 V4] = svd(h4);
h_down1(:,:) = downlink(i,:,1:8,9);
h_down2(:,:) = downlink(i,:,9:16,9);
h_down3(:,:) = downlink(i,:,17:24,9);
h_down4(:,:) = downlink(i,:,25:32,9);
%prompt1 = "What is the number of effective antenna? ";
x1 = x/4; %input(prompt1);
ef1 = 8 - x1 + 1;
p_self1 = V1(:,ef1:8); % Softnull precoding
p_self2 = V2(:,ef1:8);
p_self3 = V3(:,ef1:8);
p_self4 = V4(:,ef1:8);
Heff1 = h_down1*p_self1;
Heff2 = h_down2*p_self2;
Heff3 = h_down3*p_self3;
Heff4 = h_down4*p_self4;
%%
pdown1 = Heff1'*((Heff1*Heff1')^(-1)); 
pdown2 = Heff2'*((Heff2*Heff2')^(-1)); 
pdown3 = Heff3'*((Heff3*Heff3')^(-1)); 
pdown4 = Heff4'*((Heff4*Heff4')^(-1)); 
s_down = ones(4,1);
mp1 = (p_self1*pdown1*s_down);  
mp2 = (p_self2*pdown2*s_down);
mp3 = (p_self3*pdown3*s_down);
mp4 = (p_self4*pdown4*s_down);
x_down1 = (1/(norm(mp1,'fro')^2))*mp1;
x_down2 = (1/(norm(mp2,'fro')^2))*mp2;
x_down3 = (1/(norm(mp3,'fro')^2))*mp3;
x_down4 = (1/(norm(mp4,'fro')^2))*mp4;
SI_term1 = h1*x_down1; 
SI_term2 = h2*x_down2;
SI_term3 = h3*x_down3;
SI_term4 = h4*x_down4;
SI_power1 = norm(SI_term1,'fro')^2;
SI_power2 = norm(SI_term2,'fro')^2;
SI_power3 = norm(SI_term3,'fro')^2;
SI_power4 = norm(SI_term4,'fro')^2;
z4 = db(SI_power1+SI_power2+SI_power3+SI_power4);
sub4 = z4 - SI_power;
per4 = (-sub4*100)/SI_power;

%%

H1(:,:) = internal_ip2(i,33:64,1:16,9);
H2(:,:) = internal_ip2(i,33:64,17:32,9);
[u1 s1 v1] = svd(H1);
[u2 s2 v2] = svd(H2);
H_down1(:,:) = downlink(i,3:4,1:16,9);
H_down2(:,:) = downlink(i,3:4,17:32,9);
x2 = x/2;
ef2 = 16 - x2 + 1;
P_self1 = v1(:,ef2:16); % Softnull precoding
P_self2 = v2(:,ef2:16);
heff1 = H_down1*P_self1;
heff2 = H_down2*P_self2;
Pdown1 = heff1'*((heff1*heff1')^(-1)); 
Pdown2 = heff2'*((heff2*heff2')^(-1)); 
sdown = ones(2,1);
Mp1 = (P_self1*Pdown1*sdown);  
Mp2 = (P_self2*Pdown2*sdown);
X_down1 = (1/(norm(Mp1,'fro')^2))*Mp1;
X_down2 = (1/(norm(Mp2,'fro')^2))*Mp2;
sI_term1 = H1*X_down1; 
sI_term2 = H2*X_down2;
sI_power1 = norm(sI_term1,'fro')^2;
sI_power2 = norm(sI_term2,'fro')^2;
z2 = db(sI_power1+sI_power2);
sub2 = z2 - SI_power;
per2 = (-sub2*100)/SI_power;

%%

H(:,:)=internal_ip2(i,33:64,1:32,9);
[u s v] = svd(H);
h1_down(:,:) = downlink(i,1,1:32,9);
x3 = x;
ef3 = 32 - x3 + 1;
P1_self = v(:,ef3:32); % Softnull precoding
Heff_1 = h1_down.'*P1_self;
Pdown_1 = Heff_1'*((Heff_1*Heff_1')^(-1)); 
MP_1 = (P1_self*Pdown_1);
X_down_1 = (1/(norm(MP_1,'fro')^2))*MP_1;
SI1_term = H*X_down_1; 
SI1_power = db(norm(SI1_term,'fro')^2);
sub1 = SI1_power - SI_power;
per1 = (-sub1*100)/SI_power;

%% calculating SVD and SI
function SI = calculate_SVD(SI_matrix,M_Tx,D_Tx)
    [U,S,V] = svd(SI_matrix);
    for i=1:1:M_Tx-D_Tx
        S(i,i) = 0;
    end
    H_soft = U*S*V';
    % rescaling due to using fewer antenna
    H_soft = H_soft*(M_Tx/D_Tx);
    SI = norm(H_soft,'fro')^2;
    
    % calculating SI-Sotfnull
    %SI1 = sum(H_soft.*conj(H_soft),"all");
   % a = H_soft*H_soft';
   % SI=0;
   % for k = 1:dim
    %    l=sum(a(k,k:dim)+a(k:dim,k).');
    %    m=l-a(k,k);
     %   SI=SI+m;
   % end
end

%% calculate the capacity 
% Define the number of antennas at the base station and clients
function [capacity_uplink, capacity_downlink, capacity_all] = cap(H_uplink, H_downlink, H_SI)
    nT_BS = 28;
    nR_BS = 28;
    nT_clients = 1;
    nR_clients = 1;
    
    % Define the number of clients
    nClients = 4;
    
    
    % Define the uplink and downlink channel matrices for each client
    %H_uplink = randn(nR_BS, nT_clients, nClients);
    %H_downlink = randn(nR_clients, nT_BS, nClients);
    
    % Define the self-interference channel matrix
    %H_SI = randn(nR_BS, nT_BS);
    
    % Compute the covariance matrix of the transmitted signals and the self-interference for each client
    R = eye(nT_clients); % Assume the transmitted signals have unit covariance
    capacity_uplink = 0;
    for i = 1:nClients
    capacity_uplink = capacity_uplink + log2(det(eye(nR_BS) + 100*H_uplink(:,:,i) * H_uplink(:,:,i)'))-log2(det(eye(nR_BS)+(H_SI * H_SI')));
    end
    
    % Compute the capacity using the mutual information formula
    capacity_downlink = 0;
    for i = 1:nClients
    capacity_downlink = capacity_downlink + log2(det(eye(nT_BS) + 100* H_downlink(:,:,i)' *  H_downlink(:,:,i)));
    end
    capacity_all = capacity_uplink +  capacity_downlink;

%% Here's an example code in Matlab to compare the self-interference and capacity of a hybrid beamforming
% system with different numbers of quantization bits:
% The "quantizeCoefficients" function in the below code should be implemented by the user and 
% it should perform quantization of the beamforming coefficients with a given resolution. The code 
% below is just an example and can be adjusted according to specific requirements and conditions.

% Define the number of antennas at the base station and clients
nT_BS = 4;
nR_BS = 4;
nT_clients = 2;
nR_clients = 2;

% Define the number of clients
nClients = 3;

% Define the uplink and downlink channel matrices for each client
H_uplink = randn(nR_BS, nT_clients, nClients);
H_downlink = randn(nR_clients, nT_BS, nClients);

% Define the self-interference channel matrix
H_SI = randn(nR_BS, nT_BS);

% Define the transmitted signals covariance matrix
R = eye(nT_clients);

% Define the number of quantization bits
numBits = [2 4 6 8];

% Preallocate the self-interference and capacity arrays
selfInterference = zeros(length(numBits), 1);
capacity = zeros(length(numBits), 1);

% Loop over the different number of quantization bits
for i = 1:length(numBits)
% Define the quantization resolution
quantizationResolution = 2^(-numBits(i));

% Calculate the quantized beamforming coefficients
W_quantized = quantizeCoefficients(W, quantizationResolution);

% Compute the covariance matrix of the transmitted signals and the self-interference
P = H_uplink * R * H_uplink' + H_SI * W_quantized * W_quantized' * H_SI';

% Compute the self-interference
selfInterference(i) = trace(W_quantized * W_quantized' * H_SI' * H_SI);

% Compute the capacity
capacity(i) = 0;
for j = 1:nClients
    capacity(i) = capacity(i) + log2(det(eye(nR_clients) + H_downlink(:,:,j) * P * H_downlink(:,:,j)'));
end
end

% Plot the self-interference and capacity versus the number of quantization bits
figure;
plot(numBits, selfInterference, '-o');
xlabel('Number of quantization bits');
ylabel('Self-interference (dB)');

figure;
plot(numBits, capacity, '-o');
xlabel('Number of quantization bits');
ylabel('Capacity (bits/s/Hz)');






%% Note Here's an example code in Matlab to compare the self-interference and capacity of a hybrid beamforming
% system with different numbers of quantization bits:
% The "quantizeCoefficients" function in the below code should be implemented by the user and 
% it should perform quantization of the beamforming coefficients with a given resolution. The code 
% below is just an example and can be adjusted according to specific requirements and conditions.

% Define the number of antennas at the base station and clients
nT_BS = 4;
nR_BS = 4;
nT_clients = 2;
nR_clients = 2;

% Define the number of clients
nClients = 3;

% Define the uplink and downlink channel matrices for each client
H_uplink = randn(nR_BS, nT_clients, nClients);
H_downlink = randn(nR_clients, nT_BS, nClients);

% Define the self-interference channel matrix
H_SI = randn(nR_BS, nT_BS);

% Define the transmitted signals covariance matrix
R = eye(nT_clients);

% Define the number of quantization bits
numBits = [2 4 6 8];

% Preallocate the self-interference and capacity arrays
selfInterference = zeros(length(numBits), 1);
capacity = zeros(length(numBits), 1);

% Loop over the different number of quantization bits
for i = 1:length(numBits)
% Define the quantization resolution
quantizationResolution = 2^(-numBits(i));

% Calculate the quantized beamforming coefficients
W_quantized = quantizeCoefficients(W, quantizationResolution);

% Compute the covariance matrix of the transmitted signals and the self-interference
P = H_uplink * R * H_uplink' + H_SI * W_quantized * W_quantized' * H_SI';

% Compute the self-interference
selfInterference(i) = trace(W_quantized * W_quantized' * H_SI' * H_SI);

% Compute the capacity
capacity(i) = 0;
for j = 1:nClients
    capacity(i) = capacity(i) + log2(det(eye(nR_clients) + H_downlink(:,:,j) * P * H_downlink(:,:,j)'));
end
end

% Plot the self-interference and capacity versus the number of quantization bits
figure;
plot(numBits, selfInterference, '-o');
xlabel('Number of quantization bits');
ylabel('Self-interference (dB)');

figure;
plot(numBits, capacity, '-o');
xlabel('Number of quantization bits');
ylabel('Capacity (bits/s/Hz)');



%% This part is extra for calculating optimum weights for phase shifters.

load('trace-internal-meas-2023-1-14-12-53-49_1_30_0_29_chanA');
h(:,:) = user_csi_internalA(1,:,:,1);
N_Tx = 30;
N_Rx = 30; 
phi = 45;
phi = deg2rad(phi);
d_h = 0.07935;
d_v = 0.06668;
P = 1;
d_g = sqrt(P) * N_Tx;
freq = 3.5 * 10^9;
lambda = 3 *10^8 / freq;
w = zeros(N_Tx,1);
for aa = 1:8
    for bb = 1:4
        countup(aa,bb) = sqrt((bb*d_h)^2 + (aa*d_v)^2); 
    end
end
a=reshape(countup',1,32);
a=a(1:30);
omega = (1i*2*pi*a*cos(phi))/lambda;
a_k = exp(omega);
D = zeros(N_Rx + 1,1);
D(1,1) = d_g;
B = cat(1,a_k,h);
D_prim = cat(1,real(D),imag(D));
B_prime = cat(2,(cat(1,real(B),imag(B))),(cat(1,-(imag(B)),real(B))));
w_prime = cat(1,real(w),imag(w));
Q_prime = zeros(N_Tx,2*N_Tx,2*N_Tx);
for cc = 1:N_Tx
    for dd = 1:2*N_Tx
        if cc == dd || cc + N_Tx == dd
            Q_prime(cc,dd,dd) = 1;
        end
    end
end
w_prime2 = cat(1,w_prime,1);
B_prime2 = cat(2,(cat(1,B_prime'*B_prime,D_prim'*B_prime)),(cat(1,-(B_prime'*D_prim),(vecnorm(D_prim))^2)));
Q_prime2 = zeros(N_Tx+1,2*N_Tx,2*N_Tx);

for ee = 1:N_Tx+1
    for ff = 1:2*N_Tx+1
        if (ee ~= N_Tx+1 && (ff == ee || ff == ee + N_Tx)) || (ee == N_Tx+1 && ff == 2*N_Tx+1)
            Q_prime2(ee,ff,ff) = 1;
        end
    end
end
Q = cell(1,31);
for n=1:31
q(:,:) = Q_prime2(n,:,:);
Q{1,n}=q(:,:);
end

%% SDP optimization for W
cvx_begin 
    variable W(2*N_Tx+1,2*N_Tx+1) semidefinite 
    minimize( trace(B_prime2*W) )
    subject to
    for i = 1:30
        tr(:,:) = Q_prime2(i,:,:)
        trace(tr*W) == P;
    end
    tr1(:,:) = Q_prime2(31,:,:)
    trace(tr1*W) == 1;
    W >= 0 ;
cvx_end
%% match to angles based on num of quantization bits
di1 = diag(W);
di = sqrt(di1);
temp1 = (length(di)-1)/2;
for r = 1:temp1
    w(r) = complex(di(r),di(r+temp1)) ;
end
[theta,rho] = cart2pol(di(1:temp1),di(temp1+1:2*temp1));
numBits = [2 3 4 5 6 7 8];
for i = 1:length(numBits)
    for j = 1:2^(numBits(i))
    quantizationResolution(i,j) =2*pi*(2^(-numBits(i))*(j-1));
    end
end

for i = 1:length(numBits)
    for k = 1:temp1
        for j = 1:2^(numBits(i)) 
            diff(j) = abs(theta(k) - quantizationResolution(i,j));   
        end
        [min_diff_j, min_index_j] = min(diff);
        w_q(k) = quantizationResolution(i,min_index_j);
    end
    w_qu(:,i) = w_q;
end
for i = 1:length(numBits)
    for j = 1:temp1
    [x y] = pol2cart(w_qu(j,i),rho(j));
    w_opt(j,i) = complex(x,y);
    end
end

%% Generate a PSD matrix A of size N by N


M = W; 

% Compute the eigendecomposition of M
[Q, Lambda] = eig(M);

% Compute the vector v of square roots of eigenvalues
v = sqrt(diag(Lambda));

% Construct the product v*v^T
M_hat = v * v';

% Check that M_hat is equal to M up to numerical precision
norm(M - M_hat)
%% calcuate the weight
temp = v;
for r = 1:(length(temp)-1)/2
    z(r) = complex(temp(r),temp(r+((length(temp)-1)/2))) ;
end

%% norm nuclear rank minimization

% define the positive semidefinite matrix A
A = W;
% define the rank of low rank approximation
k = 1;
% define the optimization problem
cvx_begin
    variable X(61,61) semidefinite
    minimize (norm_nuc(X))
    subject to 
        X==semidefinite(61)
        trace(X)==k
cvx_end
% get the low rank approximation of A
X_approx = X;
% calculate between the original matrix and the approxiimation
diff = A - X_approx;
% calculate the Frobenius norm of difference
Ferobenius_norm = norm(diff, 'fro');
%% Rank Constraint

% Define the positive semi-definite matrix A
A = W;

% Define the size of the matrix
n = size(A, 1);

% Define the desired rank
k = 1;

% Define the optimization variables
cvx_begin
variable X(n, n) semidefinite
minimize(norm(A - X, 'fro'))
subject to
rank(X) <= k
cvx_end

% Get the low-rank approximation of A
X_approx = X;

% Calculate the difference between the original matrix and the approximation
diff = A - X_approx;

% Calculate the Frobenius norm of the difference
Frobenius_norm = norm(diff, 'fro');

