how to to code iterative summations in matlab

(EDIT)

hello,

im trying to replicate the graph below:

to be exact, for the time being im trying to replicate the black LoS curve, plotted through coding and simulating the equation below:

in the equation,we have K angles Q0_0 with values randomized between (0;2Pi) i refered to as Q in my code ,and K random angles Q0_1 random in value between (0;2PI) refered to as W in my code,

in the first G(.,.) , G1( Q(k), Q(i) ) where i and k are the i_th and k_th Q angles

and the second G(.,.), G2( Q(k), W(i) ) where i is the i_th angle of W, and k is the k_th Q angles

due do the random nature of generating the random angles Q and W,

beta_bar and SNR0 are given in dB we convert them to linear scale (and to avoid any errors in the case of 1/SNR0).

the original code used to plot the figure above is as follows:

%%

SNR = 0;
SNR = 10^(SNR/10);

betabar = -10;
betabar = 10^(betabar/10);

K = 1:20;
M = [10 100];

Kmax = max(K);
Mmax = max(M);
numberOfRealizations = 1000;

Q_00 = 2*pi*rand(1,Kmax,numberOfRealizations);
Q_01 = 2*pi*rand(1,Kmax,numberOfRealizations);

antennaSpacing = 1/2;
SE_MR_LoS = zeros(length(K),length(M));

%% 
for n = 1:numberOfRealizations
    disp([num2str(n) ' realizations out of ' num2str(numberOfRealizations)]);
    %Go through the range of number of UEs
    for kindex = 1:length(K)
        %Go through the range of number of BS antennas
        for mindex = 1:length(M)
            %Compute the SE with MR under LoS propagation using (1.43) for
            %one realization of the UE angles

            argumentsDesired = pi*antennaSpacing*( repmat(sin(Q_00(1,1:K(kindex),n)),[K(kindex) 1])...
                                                -  repmat(sin(Q_00(1,1:K(kindex),n))',[1 K(kindex)]) );
            argumentsInterfering = pi*antennaSpacing*( repmat(sin(Q_00(1,1:K(kindex),n)),[K(kindex) 1])...
                                                    -  repmat(sin(Q_01(1,1:K(kindex),n))',[1 K(kindex)]) );
            oneminuscos = sin(argumentsDesired).^2 + eye(K(kindex));

            %Compute the uplink SE with MR combining
            SE_MR_LoS(kindex,mindex) = SE_MR_LoS(kindex,mindex) +...
                sum(log2(1 + SNR*M(mindex)*ones(1,K(kindex)) ./...
                ( (SNR)*sum( (sin(M(mindex)*argumentsDesired)).^2 ./ (M(mindex)*oneminuscos),1) +...
                  (SNR)*betabar*sum( (sin(M(mindex)*argumentsInterfering)).^2 ./ (M(mindex)*(sin( argumentsInterfering)).^2), 1) + 1)))/numberOfRealizations;
        end
    end 
end

%% Plot the simulation results for MR combining
hold on; box on; grid on;

for mindex = 1:length(M)
    plot(K,SE_MR_LoS(:,mindex),'k-','LineWidth',1);  
end

xlabel('Number of UEs (K)');
ylabel('Average sum SE [bit/s/Hz/cell]');
ylim([0 120]);

i was told to avoid plagia i have to write my own code, so i did my best to write a code that function similarly. my best attempt code is below:

% parameters
dH = 1/2;
M = 100;
K_UE = 1:20;
Kmax = max(K_UE);
num_realizations = 1000;

SNR0_dB = 0;
SNR0 = 10^( SNR0_dB /10);

beta_bar_dB = -10;
beta_bar = 10^( beta_bar_dB /10);

% Compute SE_LOS for each value of K
TOT_SUM = 0;
SE_LOS = zeros(Kmax,1);

W = 2*pi*rand(Kmax,num_realizations);   %  W  angles of each inter_cell UE 
Q = 2*pi*rand(Kmax,num_realizations);   %  Q  angles of each intra_cell UE 

%%     
for kk = 1:length(K_UE)
    disp([num2str(kk) ' realizations out of ' num2str(Kmax)]);
    G1 = zeros(Kmax, Kmax, num_realizations); % Initialize G1 as a 3D array
    G2 = zeros(Kmax, Kmax, num_realizations); % Initialize G2 as a 3D array
    
    for j = 1:num_realizations
        
        for i = 1:kk

            for k = 1:kk
                G_intra = pi*dH*( sin( Q(i,j) ) - sin( Q(k,j)) ) ;
                G_inter = pi*dH*( sin( Q(i,j) ) - sin( W(k,j)) ) ;                          
%                 Gx(i,k,1) = 2*pi*dH*( sin( Q(i,j) ) - sin( Q(k,j)) );
                
                if  i ~= k && sin( Q(i,j)) ~= sin( Q(k,j))
                    G1(i,k,j) =  (sin( M* G_intra )) .^2./...
                                 M* (sin( G_intra )) .^2;
                    
                else if i ~= k && sin( Q(i,j)) ==  sin( Q(k,j))
                        G1(i,k,j) = M;
                    else
                        G1(i,k,j) = 0;
                    end
                end
               
                
                if  sin( Q(i,j)) ~= sin( W(k,j))                   
                    G2(i,k,j) =  (sin( M* G_inter )) .^2./...
                                 M* (sin( G_inter )) .^2;
                    
                else
                    G2(i,k,j) = M;
                end

            end   
            
        end 
        
    end
    alpha = log( 1  + M ./...  
             ((sum(G1(:)) )./num_realizations + beta_bar* (sum(G2(:)))./num_realizations + (1./SNR0))) ; 
     
    TOT_SUM = TOT_SUM +  alpha; 
    SE_LOS(kk) = TOT_SUM;   
          
end

%% Plot SE_LOS as function of K_UE
grid on; box on; hold on;
plot(K_UE,  SE_LOS(:,1));
grid on; box on;

xlabel('Number of UEs (K)')
ylabel('Sum SE (bits/s/Hz)')

i decided to take the nested loops approach, because i lack understanding when it comes to matrices,

i believe the i and k nested loops

for i = 1:kk
            for k = 1:kk
                G_intra = pi*dH*( sin( Q(i,j) ) - sin( Q(k,j)) ) ;
                G_inter = pi*dH*( sin( Q(i,j) ) - sin( W(k,j)) ) ;              

should function similarly to the repmat approach in the original code and should lead to similar results to that of the "argumentDesired" matrix considering they'll end up being summed,

            argumentsDesired = pi*dH*( repmat(sin( Q_00(1,1:K(kindex),n)),[K(kindex) 1])...
                                    -  repmat(sin( Q_00(1,1:K(kindex),n))',[1 K(kindex)]) );

before i faced issue making the SE_LoS summation add previous values to current ones, but now i don't understand why my curve doesnt resemble the desired code, as my curve is almost linear:

i also didn't understand why in the original code, why the SE_LoS is summed despite using adding previous iterations

 SE_MR_LoS(kindex,mindex) = SE_MR_LoS(kindex,mindex) + sum( full equation)
 

i hope this post wasnt too long,i apologize for the bad english, all advice and help is greatly appreciated!

[Bruno Luong]

The first code have (1-cos(...)) not in the formula nor in the second code. So clearly they don't expressthe same thing.

[Bruno Luong]

BTW in the first code

argumentsInterfering = 2*pi*antennaSpacing*( repmat(sin(varphiDesired(1,1:K(kindex),n)),[K(kindex) 1])...
    -  repmat(sin(varphiInterfering(1,1:K(kindex),n))',[1 K(kindex)]) );

is it normal that varphiDesired is mixed with varphiInterfering?

[Fares Zaritt]

hello,

regarding the use of 1-cos(...), its actually still the same equation. they for some reason decided to implement the trigonpmetry "half angle formula: sin^2(2A) = (1-cosA) /2, also multiplied the equation with ( SNR /SNR ), just to turn 1/SNR in the original formula to 1.

i rewrote a sin version of their code, and replaced all the long and fancy words with symbols to better keep track of the original formula coded:

%%

SNR = 0;
SNR = 10^(SNR/10);

betabar = -10;
betabar = 10^(betabar/10);

K = 1:20;
M = [10 100];

Kmax = max(K);
Mmax = max(M);
numberOfRealizations = 1000;

Q_00 = 2*pi*rand(1,Kmax,numberOfRealizations);
Q_01 = 2*pi*rand(1,Kmax,numberOfRealizations);

antennaSpacing = 1/2;
SE_MR_LoS = zeros(length(K),length(M));

%% 
for n = 1:numberOfRealizations
    disp([num2str(n) ' realizations out of ' num2str(numberOfRealizations)]);
    %Go through the range of number of UEs
    for kindex = 1:length(K)
        %Go through the range of number of BS antennas
        for mindex = 1:length(M)
            %Compute the SE with MR under LoS propagation using (1.43) for
            %one realization of the UE angles

            argumentsDesired = pi*antennaSpacing*( repmat(sin(Q_00(1,1:K(kindex),n)),[K(kindex) 1])...
                                                -  repmat(sin(Q_00(1,1:K(kindex),n))',[1 K(kindex)]) );
            argumentsInterfering = pi*antennaSpacing*( repmat(sin(Q_00(1,1:K(kindex),n)),[K(kindex) 1])...
                                                    -  repmat(sin(Q_01(1,1:K(kindex),n))',[1 K(kindex)]) );
            oneminuscos = sin(argumentsDesired).^2 + eye(K(kindex));

            %Compute the uplink SE with MR combining
            SE_MR_LoS(kindex,mindex) = SE_MR_LoS(kindex,mindex) +...
                sum(log2(1 + SNR*M(mindex)*ones(1,K(kindex)) ./...
                ( (SNR)*sum( (sin(M(mindex)*argumentsDesired)).^2 ./ (M(mindex)*oneminuscos),1) +...
                  (SNR)*betabar*sum( (sin(M(mindex)*argumentsInterfering)).^2 ./ (M(mindex)*(sin( argumentsInterfering)).^2), 1) + 1)))/numberOfRealizations;
        end
    end 
end

%% Plot the simulation results for MR combining
hold on; box on; grid on;

for mindex = 1:length(M)
    plot(K,SE_MR_LoS(:,mindex),'k-','LineWidth',1);  
end

xlabel('Number of UEs (K)');
ylabel('Average sum SE [bit/s/Hz/cell]');
ylim([0 120]);

ill edit my question to include this code now

also regarding the varphiDesired being mixed with varphiInterfering, now renamed to Q_00 and Q_01 respectively, the short answer is yes.

( K is the number of antennas and each antenna has it's own angle , the context for all of this is imagine two rooms next to each other, room 0 and room 1, with an internet router inside both rooms, each router serves K phones/users. now say youre using the router in room 0, all K phones in room 1 will cause interference on the routers as well as all the phones with you in room 0. and are the angles at which the router recieved the phone signal)

i apologize for not doing my due diligence explaining my question

[Bruno Luong]

Just wonder, if you have both codes why can't you just run and see what is the difference?

[Fares Zaritt]

i edited the code in my previous reply because it was faulty,

as to your last question, i was told to write a code that produces similar results to that code, but when running my code and the original code, i cant identify why my code doesnt produce similar results

if it was up to me, i'd avoid the headache and just use the original code with and add a reference to the original creater, but i was told this will still be considered plagia for some reason...