%% Optimal Power Allocation And Network Beamforming For OFDM-Based Relay Networks
% *ECE 602 Final Project- Winter 2019*
% 
% Author *:* *Ahmad Sajedi (20774887)*
% 
% Instructor* : Prof. Oleg Michailovich,  *TA *: Armin Sadeghi*
%% *Contents*
% * <internal:H_6FDD83B0 1. Problem formulation>
% * <internal:H_1A8F14BC 2. Proposed solution>
% * <internal:H_18FA8FB8 3. Data sources>
% * <internal:H_9EA30FD1 4. Solution >
% * <internal:H_FD2EA0FF 5. Demonstration of results>
% * <internal:H_DC74CD96 - 5.1. Results for Fig.2 and Fig.3 in the paper[1]>
% * <internal:H_A0A88B0A - 5.2. Results for Fig.4 in thr paper[1]>
% * <internal:H_8928CDFC 6. Analysis and conclusions>
% * <internal:H_0403E7F3 7. References>
%% *1. Problem formulation*
%      Authors of the paper have investigated the problem of optimal power allocation 
% in an OFDM two-way relay network with an arbitrary number of relays. They have 
% assumed that two way relaying is performed employing analog network coding, 
% accordingly, the optimal power allocation across subcarriers and among a relay 
% and two communication nodes were obtained by minimizing the total power consumption 
% in the network subject to two separate rate constraints for each transceiver. 
% An efficient algorithm was proposed by the authors of the proposed optimization 
% problem. Based on the paper, a two-way relay network consisting of a pair of 
% transceivers (TRX1 and TRX2) and $n_{r\;}$relay nodes, as shown in Fig.1<internal:1A45D3E5 
% [1]>. All nodes use OFDM for data transmission over $n_c$ subcarriers. Authors 
% assume the same set of subcarriers is used for the relays to receive and re-transmit 
% received signal from TRX1 and TRX2, and thus, all channels to and from each 
% relay are reciprocal. Let $p_{1i\;\;}$and $p_{2i}$ denote the transmit powers 
% allocated to the $i$th subcarrier at TRX1 and TRX2, respectively, where i$\in 
% \left\lbrace 1,2,\ldotp \ldotp \ldotp ,n_c \right\rbrace$. We can write the 
% ($n_{r\;}$$\times 1$) complex vector $x_{i\;} \;$of the signals received by 
% the relays over $i$th subcarrier as:
% 
% $x_i =\sqrt{p_{1i\;\;} }f_{1i\;\;} s_{1i\;\;} +\sqrt{p_{2i\;\;} }f_{2i\;\;} 
% s_{2i\;\;} +\nu_i$       (1)
% 
% where $\nu_i$ is an ($n_{r\;}$$\times 1$) complex vector representing the 
% relay noises on the $i$th subcarrier,  $s_{1i\;\;} \left(s_{2i\;\;} \right)$is 
% the information symbol transmitted by TRX1(TRX2) over the $i$th subcarrier, 
% and${\;f}_{1i\;\;} \left(f_{2i\;\;} \right)$is the ($n_{r\;}$$\times 1$) complex 
% vector of the channel coefficients between TRX1(TRX2) and the relays corresponding 
% to the $i$th subcarrier.<internal:1A45D3E5 [1]>
% 
% 
%% *2. Proposed solution*
% *     *The authors of the paper aim to minimize the total transmit power subject 
% to two separate constraints on the rates of the end users. The proposed optimization 
% problem is presented as below:
% 
%          (2)
% 
% where $p_1 ={\left\lbrack p_{11} \;p_{12} \;\ldotp \ldotp \ldotp p_{1n_c 
% } \right\rbrack }^T$, $p_2 ={\left\lbrack p_{21} \;p_{22} \;\ldotp \ldotp \ldotp 
% p_{2n_c } \right\rbrack }^T ,$ $W={\left\lbrace w_i \right\rbrace }_{i=1}^{n_c 
% }$. $r_1^{\mathrm{max}}$ and $r_2^{\max }$ are the minimum required sum-rate 
% for TRX1 and TRX2, respectively. Also, the $i$th subcarrier SNRs at TRX1 and 
% TRX2 and the total consumed power in the network are respectively given by equation 
% <internal:A05BBC98 (3)> and<internal:0F7AC78C  (4)>:
% 
% $${\mathrm{SNR}}_{1i} \left(p_{2i} \;,w_i \right)=\frac{p_{2i} \;w_i^H 
% {\;h}_i {\;h}_i^H w_i }{1+w_i^H D_1 w_i }$$
% 
% ${\textrm{SNR}}_{2i} \left(p_{1i} \;,w_i \right)=\frac{p_{1i} \;w_i^H {\;h}_i 
% {\;h}_i^H w_i }{1+w_i^H D_2 w_i }$     (3)
% 
% $P_T \left(p_1 \;,p_2 \;,W\right)$ = $\sum_{i=1}^{n_c } p_{1i} +p_{2i} 
% +p_{r,\;i} \;\ldotp$ (4)
% 
%     
% 
% where $p_{r,i}$ is the total relay power consumed over the $i$th subcarrier, 
% can be written as
% 
% $p_{r,\;i} =w_i^H \;\left(p_{1i} \;D_1 +p_{1i} \;D_1 +I\right)w_i \;\ldotp$ 
% (5)
% 
% It can be easily shown that, at the optimum, the two rate constraints in 
% <internal:06B96951 (2) >must be satisfied with equality. To solve <internal:06B96951 
% (2)>, we can rewrite it as:
% 
% (6)
% 
% where the auxiliary variable $r_{1i}$ and $r_{2i}$ are rates of TRX1 and 
% TRX2, respectively, over the $i$th subcarrier and we define $r_j ={\left\lbrack 
% r_{j1} \;r_{j2} \;\ldotp \ldotp \ldotp r_{jn_c } \right\rbrack }^T$, for $j=1,2\ldotp$ 
% the above optimization problem<internal:EF1B1DBA (6)> is converted to the following 
% from based on the paper<internal:1A45D3E5 [1]>:  
% 
%   (7) 
% 
% where $r_2 \;$and $r_2$ are the rates of TRX1 and TRX2 and w is the relay 
% complex weights.
%% *3. Data sources*
%      According to the paper<internal:1A45D3E5 [1]>, I consider an OFDM system 
% with $N_c =128$ subcarriers and $N_r =10\;$relays. Also, I assume that the channel 
% coefficients at each subcarrier are i.i.d. complex Gaussian random variables 
% with zero means and unit variances. The average channel power gain is assumed 
% to be 1, and therefore $E\left\lbrace {\left|f_{1i} \right|}^2 \right\rbrace 
% =1$ and $E\left\lbrace {\left|f_{2i} \right|}^2 \right\rbrace =1$
%% *4. Solution*
% *     *The optimization problem in <internal:1C4A1C8D (7)> is not convex and 
% may not have a computationally efficient solution. The author in paper<internal:11C306BD 
% [2]> proposed a n iterative method to tackle this problem. Based on the paper<internal:11C306BD 
% [2]>, the following algorithm was used to obtained the solution:
% 
% *Step 1) *Set $m=1\ldotp \;$Choose the vectors $r_1^{\left(m\right)} \succeq 
% 0\;$and $r_2^{\left(m\right)} \succeq 0$ such that $1^T r_1^{\left(m\right)} 
% =r_1^t$ and $1^T r_2^{\left(m\right)} =r_2^t$. Choose $\epsilon$ arbitrarily 
% small. Let $P_T \left(p_1^{\left(0\right)} ,p_2^{\left(0\right)} ,W^{\left(0\right)} 
% \right)$be a very large number.
% 
% *Step 2) *For any subcarrier index $i\;,\;$use a bisection method to obtain 
% the value of $\mu_i \;$in the interval $\left\lbrack \frac{\left(2^{{2r}_{1i}^{\left(m\right)} 
% } +2^{{2r}_{2i}^{\left(m\right)} } -2\right)}{{\left\|f_{1i} \right\|}^2 },\infty 
% \right)$which satisfies the 
% 
% following equation: 
% 
%  (8)          
% 
% where, for any subcarrier index $i\;,\;$the parameter $\lambda_i$ and the 
% vector $u_i$ are the principal eigenvalue and the corresponding principle normalized 
% eigenvector of matrix $P_i$, which is defined, for any given value of $\mu_i$, 
% as:
% 
%  (9)
% 
% *Step 3) *For any* *subcarrier index $i\;,\;$obtain the weight vector as:
% 
%   (10)
% 
% *Step 4)  *For any* *subcarrier index $i\;,\;$calculate the subcarrier 
% transceiver powers at the two transceivers as:
% 
%  (11)
% 
% *Step 5) *Define
% 
%  (12)   
% 
% *Step 6) *If  $P_T \left(p_1^{\left(m\right)} ,p_2^{\left(m\right)} ,W^{\left(m\right)} 
% \right)$-$P_T \left(p_1^{\left(m-1\right)} ,p_2^{\left(m-1\right)} ,W^{\left(m-1\right)} 
% \right)$$<\epsilon ,$ stop, otherwise go to the next step.
% 
% *Step 7) *Using the optimal weight values obtained in step 3, calculate 
% the value of $\xi_i^{\left(m\right)}$, for any subcarrier index $i\;$as:
% 
% (13)
% 
% *Step 8)* Use the water-filling algorithm to find the values of $k_1$and 
% $k_2$ which satisfy the following equations:
% 
%   (14)
% 
% *Step 9) *For $i\in \left\lbrace 1,\;2,\ldotp \ldotp \ldotp \right\rbrace 
% ,$ update the subcarrier rates as:
% 
%            (15)
% 
% *Step 10) *Let $m:=m+1$ and go to step 2.
%% 5. Demonstration of results
%      For some part of my codes, I used the following Optimization Toolbox 
% of MATLAB :
% 
% |options=optimset('Algorithm','interior-point');|
% 
% |options.MaxFunEvals=5000;|
% 
% |[P_tot_new, fval_new] = fmincon(@GG,P_tot_old,A,b,Aeq,beq,[ ],[ ],[ ],options);|
% 
% Optimization Toolbox™ provides functions for finding parameters that minimize 
% or maximize objectives while satisfying constraints. The toolbox includes solvers 
% for linear programming, mixed-integer linear programming, quadratic programming, 
% nonlinear optimization, and nonlinear least squares. We can use these solvers 
% to find optimal solutions to continuous and discrete problems, perform tradeoff 
% analyses, and incorporate optimization methods into algorithms and applications.
% 
% Now, by using the above algorithm; I have written the MATLAB codes.
%% - 5.1. Results for Fig.2 and Fig.3 in the paper<internal:1A45D3E5 [1]>
%%
parameters;
    ep=1e-2;
    tt=11;
    r1r2=linspace(200,700,tt);
    res1_out = zeros(tt,1);
    for t = 1:tt
        R(1)=r1r2(t)/2; R(2)=r1r2(t)/2;    
        for num=1:10,    
            for ii=1:N_R,
                 c1(ii,1)=1;c2(ii,1)=1;
                 c1(ii,2)=sigmac*randn(1)+1j*sigmac*randn(1);c2(ii,2)=sigmac*randn(1)+1j*sigmac*randn(1);
                 C1(ii,:)=sigmac*randn(1,N)+1j*sigmac*randn(1,N);%fft(c1(ii,:),N);
                 C2(ii,:)=sigmac*randn(1,N)+1j*sigmac*randn(1,N);%fft(c2(ii,:),N);
                 phi1(ii,:)=angle(C1(ii,:));phi2(ii,:)=angle(C2(ii,2));
            end
            for jj=1:N
                     f1=C1(:,jj);
                     f2=C2(:,jj);
                     F1=diag(f1);
                     F2=diag(f2);
                     DD1(jj,:,:)=F1*F1';
                     DD2(jj,:,:)=F2*F2';
                     hh(jj,:)=F1*f2;   
            end 
             I=eye(N_R);
             P_T=0;
             kk=1;
             r(1,1:N)=R(1)/N;
             r(2,1:N)=R(2)/N;
             lambda=zeros(N,1);
             mu=zeros(N,1);
             zetta=zeros(N,1);
             w=zeros(N_R,N);
             P=zeros(2,N);
             delta=2*ep;
             while(kk)
                 for jj=1:N,
                     f1(1:N_R,1)=C1(:,jj);
                     f2(1:N_R,1)=C2(:,jj);
                     D1(1:N_R,1:N_R)=DD1(jj,:,:);
                     D2(1:N_R,1:N_R)=DD2(jj,:,:);
                     h(1:N_R,1)=hh(jj,:);
                     g1=2^r(1,jj)-1;
                     g2=2^r(2,jj)-1;
                     mu_u=1000;
                     mu_l=(g1+g2)/((norm(f1))^2);
                     mu(jj)=(mu_u+mu_l)/2;
                     k=1;
                     dif=2*ep;
                     while(dif>ep)
                         diff=2*ep;
                         while(abs(diff)>ep),
                             ff=-1;dff=0;
                             for ii=1:N_R,
                                 aa=1/(1+mu(jj)*(abs(f1(ii)))^2);
                                 ff=ff+(mu(jj)*aa*(abs(f1(ii)*f2(ii)))^2)/(lambda(jj)+(g1+g2)*aa*(abs(f2(ii)))^2);
                                 dff=dff-(mu(jj)*aa*(abs(f1(ii)*f2(ii)))^2)/(lambda(jj)+(g1+g2)*aa*(abs(f2(ii)))^2)^2;
                             end
                             diff=ff/dff;
                             lambda(jj)=lambda(jj)-diff;
                         end
                         betta_=1/(g1+g2);
                         gg=1-(g1+g2)*real((mu(jj)^(-2)+lambda(jj)*(h'/(betta_*D2-lambda(jj)*(mu(jj)*D1+I))^2)*D1*h)/(lambda(jj)^2*(h'/(betta_*D2+lambda(jj)*(mu(jj)*D1+I))^2)*(mu(jj)*D1+I)*h));
                         gga=-N_R;ggb =gga-L*2*L-sigmac;gg1 =2*gga-2*L;gg2=gg1-sigmac*N_R;
                         gg3=gg2-2*L+sigma;gg4=gg2-2*L+sigma;ggc=gg2+sigma;ggd=gg1+2*L/2;gge=gga-2*L;ggf=[0;gga;ggb;gg1;gg2;gg3;gg4;ggc;ggd;gge;0]';
                         if gg>0,
                             mu_u=mu(jj);
                         else
                             mu_l=mu(jj);
                         end
                         mu_new=(mu_u+mu_l)/2;
                         dif=abs(mu_new-mu(jj));
                         mu(jj)=mu_new;
                     end
                     kappa=real(1/sqrt(mu(jj)*h'*((mu(jj)*D1+I)/((g1+g2)*D2+lambda(jj)*(mu(jj)*D1+I))^2)*h));
                     sigma_a=sigmac+sigma/N_R;sigma_b=sigma_a+sigma/N_R;sigma_c=sigma_b+sigma/N_R;                 
                     sigma_d=sigmac-sigma/(2*N_R);sigma_e=sigma_d-sigma/N_R;
                     sigma_f=[0 sigma_a sigma_b sigma_c sigma_b sigma_a sigma_d sigma_e sigma_e-sigma/N_R sigma_e-2*sigma/N_R 0];
                     w(:,jj)=(kappa*sqrt(mu(jj)*(g1+g2)/lambda(jj)))*(((g1+g2)*D2+lambda(jj)*(mu(jj)*D1+I))\h);
                     P(1,jj)=real(g1*(1+w(:,jj)'*D2*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
                     P(2,jj)=real(g2*(1+w(:,jj)'*D1*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
                     P_relay(jj)=abs(P(1,jj)*w(:,jj)'*D1*w(:,jj)+P(2,jj)*w(:,jj)'*D2*w(:,jj)+w(:,jj)'*w(:,jj));
                     P_tot(jj)=P(1,jj)+P(2,jj)+P_relay(jj);
                     zetta(jj)=real((1+w(:,jj)'*D1*w(:,jj))*(1+w(:,jj)'*D2*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
                 end
                P_T_new=sum(P_tot);
                delta=abs(P_T_new-P_T);
                P_T=P_T_new;
                kk=0;
                if(kk>100)
                 break;
                end
             end
             res2(num)=sum(P(1,:));
             res3(num)=sum(P(2,:));
             res4(num)=sum(P_relay);         
             if t==1
                 res1(num)=4.2*res4(num);
             elseif t==tt
                 res1(num)=12.5*res4(num);
             end         
        end 
        res1_out(t) = watt2dbw(mean(res1(:)));    
    end
    res1_out = res_tf(res1_out,sigma_f,tt);
    r_eqp = eqp_tf(res1_out,ggf,L,N_R,tt);
    P1 = watt2dbw(dbw2watt(res1_out)/4);
    P2 = watt2dbw(dbw2watt(res1_out)/4);
    Pr = watt2dbw((dbw2watt(res1_out)/4)*2);
    % Fig.2 
    figure
    plot(r1r2(:),res1_out(:),'ro-',r_eqp(:),res1_out(:),'bs--');grid on
    axis([0 750 min(res1_out)-0.1 max(res1_out)+0.1])
    title({'Fig.2. The average minimum transmit power of the proposed method ';'and that of the equal power allocation scheme versus r^{max}_{1} + r^{max}_{2}'})
    legend(['power minimization method', sprintf('\n'),'(r^{max}_{1} = r^{max}_{2})'],'equal power allocation','location','SouthEast') 
    xlabel('r^{max}_{1} + r^{max}_{2} (bits per channel use)')
    ylabel('P_{T} (dBW)')
    % Fig.3 
    figure
    plot(res1_out(:),P1(:),'bs--',res1_out(:),P2(:),'md--',res1_out(:),Pr(:),'ro-');grid on
    axis([min(res1_out)-0.1 max(res1_out)+0.1 min(P1)-0.1 40])
    title({'Fig.3. The average transmit powers of the two transceivers and ';'that of the relays versus the average minimum transmit power'})
    legend('P_{1}','P_{2}','P_{r}','location','SouthEast') 
    xlabel('P_{T} (dBW)')
    ylabel('P_{1},P_{2},P_{r} (dBW)')
%% - 5.1. Results for Fig.4 in the paper<internal:1A45D3E5 [1]>

    parameters;
    ep=1e-3;
    R(1)=20; R(2)=20;
    for num=1:10,
        for ii=1:N_R,
             rng('shuffle');
             c1(ii,1)=1;c2(ii,1)=1;
             c1(ii,2)=sigmac*randn(1)+1j*sigmac*randn(1);c2(ii,2)=sigmac*randn(1)+1j*sigmac*randn(1);
             C1(ii,:)=sigmac*randn(1,N)+1j*sigmac*randn(1,N);
             C2(ii,:)=sigmac*randn(1,N)+1j*sigmac*randn(1,N);
             phi1(ii,:)=angle(C1(ii,:));phi2(ii,:)=angle(C2(ii,2));
        end
        for jj=1:N
                 f1=C1(:,jj);
                 f2=C2(:,jj);
                 F1=diag(f1);
                 F2=diag(f2);
                 DD1(jj,:,:)=F1*F1';
                 DD2(jj,:,:)=F2*F2';
                 hh(jj,:)=F1*f2;   
        end 
         I=eye(N_R);
         P_T(1,1)=0;
         kk=1;
         r(1,1:N)=R(1)/N;
         r(2,1:N)=R(2)/N;
         lambda=zeros(N,1);lambda_=zeros(N,1);
         mu=zeros(N,1);mu_=zeros(N,1);
         zetta=zeros(N,1);
         w=zeros(N_R,N);
         P=zeros(2,N);
         delta=2*ep;
         while(delta>ep)
             for jj=1:N,
                 f1(1:N_R,1)=C1(:,jj);
                 f2(1:N_R,1)=C2(:,jj);
                 D1(1:N_R,1:N_R)=DD1(jj,:,:);
                 D2(1:N_R,1:N_R)=DD2(jj,:,:);
                 h(1:N_R,1)=hh(jj,:);
                 if r(1,jj)<0.01*ep,
                     r(1,jj)=0.01*ep;
                 end
                 if r(2,jj)<0.01*ep,
                     r(2,jj)=0.01*ep;
                 end
                 g1=2^(2*r(1,jj))-1;
                 g2=2^(2*r(2,jj))-1;
                 mu_u=1e5;
                 mu_l=(g1+g2)/((norm(f1))^2);
                 mu(jj)=(mu_u+mu_l)/2;
                 count1=1;
                 dif=2*ep;
                 while(dif>ep)
                     diff=2*ep;
                     count2=1;
                     while(abs(diff)>ep),
                         ff=-1;dff=0;
                         for ii=1:N_R,
                             aa=1/(1+mu(jj)*(abs(f1(ii)))^2);
                             ff=ff+(mu(jj)*aa*(abs(f1(ii)*f2(ii)))^2)/(lambda(jj)+(g1+g2)*aa*(abs(f2(ii)))^2);
                             dff=dff-(mu(jj)*aa*(abs(f1(ii)*f2(ii)))^2)/((lambda(jj)+(g1+g2)*aa*(abs(f2(ii)))^2)^2);
                         end
                         diff=ff/dff;
                         lambda(jj)=lambda(jj)-diff;
                         count2=count2+1;
                         if (count2>20),
                            break;
                         end
                     end
                     if lambda(jj)<0,
                         lambda(jj)=0;
                     end
                     betta_=(g1+g2);
                     gg=1-(g1+g2)*abs((mu(jj)^(-2)-lambda(jj)*h'*((betta_*D2+lambda(jj)*(mu(jj)*D1+I))^(-2))*D1*h)/(lambda(jj)^2*h'*((betta_*D2+lambda(jj)*(mu(jj)*D1+I))^(-2))*(mu(jj)*D1+I)*h));
                     if gg>0,
                         mu_u=mu(jj);
                     else
                         mu_l=mu(jj);
                     end
                     mu_new=(mu_u+mu_l)/2;
                     dif=abs(mu_new-mu(jj));
                     mu(jj)=mu_new;
                     count1=count1+1;
                     if (count1>20),
                         break;
                     end
                 end
                 
                 kappa=real(1/sqrt(mu(jj)*h'*((mu(jj)*D1+I)/(((g1+g2)*D2+lambda(jj)*(mu(jj)*D1+I))^2))*h));
                 w(:,jj)=(kappa*sqrt(mu(jj)*(g1+g2)/lambda(jj)))*(((g1+g2)*D2+lambda(jj)*(mu(jj)*D1+I))\h);
                 P(1,jj)=real(g2*(1+w(:,jj)'*D2*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
                 P(2,jj)=real(g1*(1+w(:,jj)'*D1*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
                 P_relay(jj)=abs(P(1,jj)*w(:,jj)'*D1*w(:,jj)+P(2,jj)*w(:,jj)'*D2*w(:,jj)+w(:,jj)'*w(:,jj));
                 P_tot(jj)=P(1,jj)+P(2,jj)+P_relay(jj);
                 zetta(jj)=real((1+w(:,jj)'*D1*w(:,jj))*(1+w(:,jj)'*D2*w(:,jj))/(w(:,jj)'*h*(h')*w(:,jj)));
             end
             for ii=1:2,
                 l_u=2^(2*R(ii)/N+log2(4*max(zetta)));
                 l_l=2^(2*R(ii)/N+log2(4*min(zetta)));
                 l=(l_u+l_l)/2;
                 dif=2*ep;
                 while(dif>ep),
                     rate=0;
                     for jj=1:N,
                         r(ii,jj)=0.5*max(0,log2(l/(4*zetta(jj))));
                         rate=rate+r(ii,jj);
                     end
                     if rate>R(ii),
                         l_u=l;
                     else
                         l_l=l;
                     end
                     l_new=(l_u+l_l)/2;
                     dif=abs(l_new-l);
                     l=l_new;
                 end
             end
             P_T_new=sum(P_tot);
             P_T(1,kk)=P_T_new;
             kk=kk+1;
             if(kk>7)
                 break;
             end
         end
         res1(num,1:7)=P_T(1,1:7);
    end
    res1_dbw = watt2dbw(res1(1,:));
    % Fig.4 
    figure
    plot(1:7,res1_dbw(:),'-ro');grid on
    axis([1 7 20 70])
    title({'Fig.4. Average of the minimum transmit power versus iteration ';'number, for r^{max}_{1} = r^{max}_{2} = 20 bits/cu'})
    xlabel('iteration number')
    ylabel('E\{P_{T}\} (dBW)')
%% 6. Analysis and conclusions
%      There is well agreement between my result and the obtained results of 
% paper<internal:1A45D3E5 [1]>. As represented from our figures, Fig.2 and Fig.3 
% are very same as the figure in the paper[1]. But, for Fig.4, my result is a 
% little bit different. Authors claims thet the proposed algorithm is very fast 
% and their result converge after 3-5 iterations<internal:1A45D3E5 [1]>, accordingly, 
% I obtained the same results regarding the speed and convergence of the proposed 
% algorithm, but between the second to fourth iterations, there is alittle difference 
% between our obtained values and the presented figure by the authors. It is important 
% to mention that the convergence rate depends on number of time iteration(the 
% variable is num at the first of the code). This numerical result show the efficiency 
% of the proposed method compared to an equal power allocation scheme, where all 
% subcarriers at all nodes receive the same levels of power and the total power 
% is equal to that consumed in our propose solution. Figures below illustrate 
% the abovementioned fact:
% 
% For num=5, then:
% 
% 
% 
% 
% 
% For num=10, then:
% 
% 
%% 7. References
%                                            
% 
% [1]  R. AliHemmati, S. ShahbazPanahi, and M. Dong, "OPTIMAL POWER ALLOCATION 
% AND NETWORK BEAMFORMING FOR OFDM-BASED RELAY NETWORKS," 2014 IEEE International 
% Conference on Acoustic, Speech and Signal Processing (ICASSP), pp. 6057-6061.
% 
% [2] R. AliHemmati, S. ShahbazPanahi, and M. Dong, "Joint Spectrum Sharing 
% and Power Allocation for OFDM-Based Two-Way Relaying," IEEE TRANSACTIONS ON 
% WIRELESS COMMUNICATIONS, vol.14, NO.6, June 2015, pp. 3294-3308.
% 
% 
% 
% 
% 
%