IA2/calculate_beamformers2.cpp

// Copyright 2011-2013, Per Zetterberg, KTH Royal Institute of Technology
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//

#include "IA2_node.hpp"
#include <itpp/itcomm.h>
#include <itpp/stat/misc_stat.h>

void IA2_node::calculate_beamformers(void) {
 
  // Perform the distrinuted IA algo found in "Approaching the Capacity
  // of Wireless Networks through Distributed Interference
  // Alignment", by Krishna Gomadam, Viveck R. Cadambe and
  // Syed A. Jafar. IEEE Globecom 2008 and ArXiv http://arxiv.org/abs/0803.3816.
  // The ArXiv version is most detailed.

  // If max_SINR=true (see below) then  the Max-SNIR version is used.

  /* The matrix Vbig[i1](:,:) contains, as it's i2:th column, the
    precoder for stream i1 on the ix_all(i2) subcarrier.
    The vector ix_all contains all the subcarrier which
    are used (not all are). 
    The subcarrier numbers we use are actually
    1,...,ai[0].rightmost_carrier
    and ai[0].leftmost_carrier to Nfft-1. These are the ones
    closest to the 0Hz unused DC subcarrier.
    Note that subcarrier numbers greater than Nfft/2 wraps
    into the negative part of the spectrum.
    The two carriers ix_all(ai[0].pilot_carrier1_sub) and 
    ix_all(ai[0].pilot_carrier2_sub)
    are in the middle of the positive and negative part of the
    spectrum, respectively. The names indicate that there would be only
    pilots being transmitted on these subcarriers. However, this is
    not the case. Since we have set "do_use_pilot_carriers=false".
    
    Note that the paper by Gadambe&Jafar referenced above
    is defined for a single carrier. We use the fact that we
    have multiple subcarriers to speed up convergence,
    The idea is that adjacent subcarriers should have similiar
    channel responses and thus we use the adjacent subcarrier
    as the starting point for the iteration. 
  */
  
  /* Define a matrix to store all the beamformers on a singe subcarrier */
  cmat v_old(num_ant_bs,num_streams);
  
  /* Initialize v_old at the beamformers we used in previous time-slot
     on the carrier pilot_carrier1_sub. */
  for (uint32_t i1=0;i1<num_streams;i1++) {
    for (uint32_t i3=0;i3<num_ant_bs;i3++) {
        v_old(i3,i1)=Vbig[i1](i3,ai[0].pilot_carrier1_sub);
    };
  };

  /* Update the beamformers starting from subcarrier index 
     ix_all(ai[0].pilot_carrier1_sub) and ending at 
     ix_all(ai[0].rightmost_carrier_sub) */
  iteration(v_old,ai[0].pilot_carrier1_sub,ai[0].rightmost_carrier_sub,1);
  /* Update the beamformers starting from subcarrier index 
     ix_all(ai[0].pilot_carrier1_sub) and ending at ix_all(0) */
  iteration(v_old,ai[0].pilot_carrier1_sub,0,-1);


  /* Initialize v_old at the beamformers we used in previous time-slot
     on the carrier pilot_carrier1_sub. */
  for (uint32_t i1=0;i1<num_streams;i1++) {
    for (uint32_t i3=0;i3<num_ant_bs;i3++) {
        v_old(i3,i1)=Vbig[i1](i3,ai[0].pilot_carrier2_sub);
    };
  };

  /* Update the beamformers starting from subcarrier index 
     ix_all(ai[0].pilot_carrier2_sub) and ending at ix_all(length_ix_all-1) 
     i.e. the maximum used subcarrier number.*/
  iteration(v_old,ai[0].pilot_carrier2_sub,ai[0].length_ix_all-1,1);

  /* Update the beamformers starting from subcarrier index 
     ix_all(ai[0].pilot_carrier2_sub) and ending at ix_all(length_ix_all-1) 
     i.e. the maximum used subcarrier number.*/
  iteration(v_old,ai[0].pilot_carrier2_sub,ai[0].leftmost_carrier_sub,-1);      



};

void IA2_node::iteration(const cmat v_old_in,int subcarrier_ix_start,int subcarrier_ix_stop,int increment) {

  cmat B(num_ant_ms,num_ant_ms); // Interference plus noise covariance matrix   
  cmat Hi(num_ant_ms,num_ant_bs); // Intermediate channel matrix
  cmat Br(num_ant_bs,num_ant_bs); // Reverse interference plus noise matrix
  cmat v_old;

  cvec v(num_ant_bs); // Intermediate precoding vector
  cvec u(num_ant_ms); // Intermediate receiving vector
  cvec h(num_ant_ms); // Intermediate kind of channel vector
  cvec hr(num_ant_bs); // Intermediate kind of channel vector
  complex<double> noise_variance;
  cmat u_old(num_ant_ms,num_ms);
  noise_variance=5901.6;
  ivec bs_antenna_ix_for_ms_ix[3];
  bool max_SINR=true;
  cmat E, Er;
  vec d, dr;
  vec SINR(3);
  vec SINR_min="1e10 1e10 1e10";

  v_old=v_old_in;

  if (!max_SINR) {
    E.set_size(num_ant_ms,num_ant_ms);
    Er.set_size(num_ant_bs,num_ant_bs);
    d.set_length(num_ant_ms);
    dr.set_length(num_ant_bs);
  };
    
  if (scenario==1) { // IA
    // The indecies of the base-station antennas, numbered from 0 to 5,
    // used by MS 0,1 and 2 in the IA approach.
    bs_antenna_ix_for_ms_ix[0]="0 1";
    bs_antenna_ix_for_ms_ix[1]="2 3";
    bs_antenna_ix_for_ms_ix[2]="4 5";
  };
  if (scenario==2) { // COMP
    // The indecies of the base-station antennas, numbered from 0 to 5,
    // used by MS 0,1 and 2 in the CoMP approach.
    bs_antenna_ix_for_ms_ix[0]="0 1 2 3 4 5";
    bs_antenna_ix_for_ms_ix[1]="0 1 2 3 4 5";
    bs_antenna_ix_for_ms_ix[2]="0 1 2 3 4 5";
  }

  for (int ix=subcarrier_ix_start;
       ix!=subcarrier_ix_stop+increment;ix+=increment) {

    for (int iter=0;iter<10;iter++) {

      /* Calculate all receive combiners */
      for (uint32_t i1=0;i1<num_streams;i1++) {

        /* B is going to be the noise+interference covariance matrix 
           for MS i1. Start with the noise. */
        eye(num_ant_ms,B);
        B=B*noise_variance;

        for (uint32_t i2=0;i2<num_streams;i2++) {
          

          //cout << "================= \n";
          //cout << "v_old=" << v_old << endl;

          /* v is the pre-coder of stream i2 */
          v=v_old.get_col(i2);
          if (!(i1==i2)) {

            
            // Fill Hi with the MIMO channel from interfering stream i2
            // from the viewpoint of stream (i.e. MS)  i1.    

            for (uint32_t i3=0;i3<num_ant_bs;i3++) {
              for (uint32_t i4=0;i4<num_ant_ms;i4++) {

                Hi(i4,i3)=H[stream_ix_to_ms_ix(i1)][bs_antenna_ix_for_ms_ix[i2](i3)](i4,ix);
              };
            }; 

            // h is the vector channel given that stream i2 use pre-coder v.
            h=Hi*v;
            
            /* B is the total downlink interference covariance matrix for i1*/
            B=B+outer_product(h,h,true);
          }; 
          //cout << "B=" << B << ";" << endl;
       
        };

        if (max_SINR) {
        
          /* v is now the pre-coder of the desired stream */
          v=v_old.get_col(i1);
          /* Fill Hi with MIMO channel of the desired stream */
          for (uint32_t i3=0;i3<num_ant_bs;i3++) {
            for (uint32_t i4=0;i4<num_ant_ms;i4++) {
              Hi(i4,i3)=H[stream_ix_to_ms_ix(i1)][bs_antenna_ix_for_ms_ix[i1](i3)](i4,ix);
            };
          }; 
          /* Calculate MSSE */
          h=Hi*v;
          u=inv(B)*h;
          u=u/norm(u,2);
          u_old.set_col(i1,u);

          //cout << "Hd=" << Hi << ";" << endl;
          //cout << "u=" << u << ";" << endl;
 
        } 
        else 
        {
          eig_sym(B,d,E);
          u=E.get_col(0);
          u_old.set_col(i1,u); 
        };
        /* Calculate the SINR given everything is perfect */
        SINR(i1)=abs(dot(conj(u),h));
        SINR(i1)=SINR(i1)*SINR(i1)/abs((dot(conj(u),B*u)));
      };
      
 
 
      /* Calcute new downlink beamformers */
      for (uint32_t i1=0;i1<num_streams;i1++) {         

        /* Br is the reciprocal interference plus noise covariance matrix */
        eye(num_ant_bs,Br);
        Br=Br*noise_variance;
        for (uint32_t i2=0;i2<num_streams;i2++) {
          if (i1!=i2) {
            /* u is receiver filter of stream i2 - now acting as transmitter */
            u=u_old.get_col(i2);
            /* Hi is the reciprocal MIMO channel from MS i2 to BS i1 */
            for (uint32_t i3=0;i3<num_ant_bs;i3++) {
              for (uint32_t i4=0;i4<num_ant_ms;i4++) {
                Hi(i4,i3)=H[stream_ix_to_ms_ix(i2)][bs_antenna_ix_for_ms_ix[i1](i3)](i4,ix);
              };
            };
            
            /* Reciprocal vector channel with u as transmitter */
            hr=hermitian_transpose(Hi)*u; 
            Br=Br+outer_product(hr,hr,true);
          };
        };
        if (max_SINR) {
          /* u is the reciprocal beamformer of stream i1 */
          u=u_old.get_col(i1);
          /* Hi is the desired reciprocal channel of stream i1 */
          for (uint32_t i3=0;i3<num_ant_bs;i3++) {
            for (uint32_t i4=0;i4<num_ant_ms;i4++) {
              Hi(i4,i3)=H[stream_ix_to_ms_ix(i1)][bs_antenna_ix_for_ms_ix[i1](i3)](i4,ix);
            };
          };

          /* Calculate new reciprocal beamformer */ 
          hr=hermitian_transpose(Hi)*u;
          v=inv(Br)*hr;
          v=v/norm(v,2);
          Vbig[i1].set_col(ix,v);
          v_old.set_col(i1,v);

          //cout << "Hd=" << Hi << ";" << endl;
          //cout << "v=" << v_old << ";" << endl;
        }
        else 
        {
          eig_sym(Br,dr,Er);
          v=Er.get_col(0);
          Vbig[i1].set_col(ix,v);
          v_old.set_col(i1,v);
        };
      };
    };
    for (uint32_t i1=0;i1<num_ms;i1++) {
      // Stor the worst SNIR of stream 0,1,2 in SINR_min
      // (worst over subcarriers). 
      if (SINR(i1)<SINR_min(i1)) 
        SINR_min(i1)=SINR(i1);  
    };
  };
  //cout << "SINR_min=" << SINR_min << endl;
};