/*************************************************************************** *   chest_rcv.c  - MIMO channel estimation module  *   begin          :  03/08/05 *   Author         : Selvavinayagam Gunabalan *   emails         : selvan@kth.se ***************************************************************************//*************************************************************************** *                                Changes *                                ------- * 03/08/05 - Selvan - Begin * 03/08/14 - Selvan - Added variance of noise estimation * 04/01/12 - ineiti - added a matched filter as option * 04/03/06 - ineiti - adjusted description * 04/03/27 - ineiti - added a very simple mimo-reception * 04/04/08 - ineiti - added a complex stats-var for testing *//*************************************************************************** *                                                                         * *   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 2 of the License, or     * *   (at your option) any later version.                                   * *                                                                         * ***************************************************************************//** *  Using the midamble the channel is estimated for the given antenna. This module * also removes the midamble part of the data stream. *  It does assume that there are a maximum of training_seq[].M senders, and that  * each sender has an increasing index. This means that a sender has at most * training_seq[].L channel taps. *  Each output corresponds to one sender. So, if you want to receive sender #2, * connect to output #2 */#include "spc.h"#include "sequences.h"#include "std.h"#include <math.h>#define DBG_LVL 0#define MAX_CHANNELS 4typedef struct {  //0-> MCCDMA, 1 -> Fra  int type; // 0  //Necessary for noise variance estimation, represnts no. Tx antennas  int num_of_antennas; // 1  // How many taps shall be used in the channel equalisation. If you  // want to use the signal for some decoding, this is best left to  // 0 (which means the signal won't be altered and you only get a channel  // estimation). If there is some hard-decision, put it to 1 or more  int calc_taps; // 4  // 0-> don't move the signal  // else align output to the strongest tap of channel "align"  // if calc_taps > 0, then align means that the number of taps calculated  // are 'around' the peak of the channel impulse response.  int align; // 1  // The length of the circular extension  int circ_ext; // 8}config_t;typedef struct {  // The channel of all N received antennas  block_t channel;  // The inversed channel  block_t channel_inv;  // The midamble  block_t midamble;  // The SNR in db  double snr;  // The general noise variance = Realpart+Imagpart  int noise_var;  //Real part  int noise_var_real;  //Imag part  int noise_var_imag;  // The length of the channel impulse response  int ch_length;  // The peak of the midamble  int peak_pos;  // peak-amplitude of the midamble  int peak_amp;  // mean midamble-amplitude  int mid_amp;  // Test  complex double var11;}stats_t;typedef struct {  SYMBOL_COMPLEX channel[ MAX_NO_OF_ANT * MAX_NO_OF_TAPS ];  SYMBOL_COMPLEX channel_inv[ MAX_NO_OF_ANT * MAX_NO_OF_TAPS ];  SYMBOL_COMPLEX *mid;  int calc_taps;  int align;  int index;}private_t;/* * The initialisation function, or constructor,  * is called the first time this module is instantiated. */int rcv_mimo_init( swr_sdb_t *context ) {  // Begin system-definitions {  config_t *config;  stats_t *stats;  MOD_INC_USE_COUNT;  if ( sizeof( private_t ) > 0 )    context->private_data = swr_malloc( sizeof( private_t ) );  swr_sdb_get_config_struct( context->id, (void**)&config );  swr_sdb_get_stats_struct( context->id, (void**)&stats );  // } End of system-definitions  config->type = 0;  config->num_of_antennas = 1;   //SISO case!  config->calc_taps = 4;  config->align = 1;  config->circ_ext = 8;  private->mid = NULL;  stats->channel.data = private->channel;  stats->channel.size = MAX_NO_OF_ANT * MAX_NO_OF_TAPS;  stats->channel.type = SIG_SYMBOL_COMPLEX;  stats->channel_inv.data = private->channel_inv;  stats->channel_inv.size = MAX_NO_OF_ANT * MAX_NO_OF_TAPS;  stats->channel_inv.type = SIG_SYMBOL_COMPLEX;  stats->midamble.size = 0;  stats->midamble.type = SIG_SYMBOL_COMPLEX;  stats->ch_length = 0;  stats->snr       = 1.0;    stats->noise_var    = 0;  stats->noise_var_real = 0;  stats->noise_var_imag = 0;  stats->peak_pos = 0;  stats->mid_amp = 0;  stats->peak_amp = 0;  stats->var11 = 0;  // Begin system-definitions  swr_sdb_free_stats_struct( context->id, (void**)&stats );  swr_sdb_free_config_struct( context->id, (void**)&config );  return 0;  // End system-definitions}/* * Every time modules from the outside change the value of a configuration parameter, * this function is called. */int rcv_mimo_reconfig( swr_sdb_t *context ) {  // Definition of variables - don't touch  config_t *config;  stats_t *stats;  swr_sdb_get_config_struct( context->id, (void**)&config );  swr_sdb_get_stats_struct( context->id, (void**)&stats );  if ( ( config->type >= TRAINING_SEQUENCES ) ||       ( config->type < 0 ) ) {    PR_DBG( 0, "Midamble-type %i not available, "	    "only 0..%i are valuable types.\n",            config->type, TRAINING_SEQUENCES - 1 );    config->type = TRAINING_SEQUENCES - 1;  }  private->calc_taps =     min( config->calc_taps, training_seq[ config->type ].L );  private->align = config->align;  if ( private->mid ){    swr_free( private->mid );    private->mid = 0;  }  stats->midamble.size = training_seq[ config->type ].Nt;  private->mid = swr_malloc( sizeof( SYMBOL_COMPLEX ) * 			     stats->midamble.size );  stats->midamble.data = private->mid;  // Definition - don't touch  swr_sdb_free_stats_struct( context->id, (void**)&stats );  swr_sdb_free_config_struct( context->id, (void**)&config );  return 0;}/* * To configure the inputs * this is called when the output-sizes change. */int rcv_mimo_configure_inputs( swr_sdb_t *context ) {  // Definition of variables - don't touch  config_t *config;  swr_sdb_get_config_struct( context->id, (void**)&config );  size_in(0) = size_out(0) + training_seq [ config->type ].Nt +    config->circ_ext;  PR_DBG( 0, "Trying to resize input, this is not good!\n" );  // Definition - don't touch  swr_sdb_free_config_struct( context->id, (void**)&config );  return 0;}/* * To configure the outputs * this is called when the input-sizes change */int rcv_mimo_configure_outputs( swr_sdb_t *context ) {  // Definition of variables - don't touch  config_t *config;  int ch, size;  swr_sdb_get_config_struct( context->id, (void**)&config );  size =  size_in(0) - training_seq[ config->type ].Nt - config->circ_ext;  for ( ch=0; ch<MAX_CHANNELS; ch++ ){    size_out(ch) = size;  }  // Definition - don't touch  swr_sdb_free_config_struct( context->id, (void**)&config );  return 0;}void do_mafi( SYMBOL_COMPLEX *in, SYMBOL_COMPLEX *out,	      complex double *f, double f_abs, int taps, int length );/* * This is the function that implements the `main method' of the class * Every class has got just ONE method/working-mode. */int rcv_mimo_pdata( swr_sdb_t *context ) {  // Definition of variables - don't touch  stats_t *stats;  config_t *config;  SYMBOL_COMPLEX *in, *out[MAX_CHANNELS], *mid, *channel, *channel_inv;  int data_len[2];  int type, num_of_antennas, peak_pos = 0, peak_amp, mid_amp, circ_ext;  short s_column; // Number of columns of shsish (S Hermitian S )inverse S hermitian matrix  short s_row;  //Number of rows of shsish (S Hermitian S )inverse S hermitian matrix  int i,j, ch; //counter  double *par; //Parity matrix = shsish  short int *seq;  short channel_offset;  short seq_offset;  short length;  double double_real = 0.;  double double_imag = 0.;  double tmp =0.;    int noise_term_real  = 0;  int noise_term_imag  = 0;  int real_term = 0;  int imag_term = 0;  double var_noise_term_imag = 0,    var_noise_term_real = 0;    int noise_var = 0;  int noise_var_real = 0;  int noise_var_imag = 0;  int seq_length, seq_ch_taps, seq_max_ant;   double sig_power = 0.;  complex double f_invs[ MAX_CHANNELS ][ 128 ];  int f_len, f_max = 0;  double f_abs = 0;  in = buffer_in(0);  for ( ch=0; ch<MAX_CHANNELS; ch++ ){    if ( port_out(ch).sdb_id >= 0 ){      PR_DBG( 4, "Got output %i\n", ch );      out[ch] = buffer_out(ch);    } else {      out[ch] = NULL;    }  }  swr_sdb_get_config_struct( context->id, (void**)&config );  type        = config->type;  seq_length  = training_seq[ type ].Nt;  seq_ch_taps = training_seq[ type ].L;  seq_max_ant = training_seq[ type ].M;  i = size_in(0) - seq_length - config->circ_ext;  // These two might actually differ...  data_len[0] = i / 2;  data_len[1] = i - data_len[0];  num_of_antennas = config->num_of_antennas;  circ_ext = config->circ_ext;  swr_sdb_free_config_struct( context->id, (void**)&config );  //Points to the last Midamble rcvd  mid = in + data_len[0] + seq_length - 1;  channel = private->channel;  channel_inv = private->channel_inv;  //Matrix Multiplication of sshish_ * received midamble  s_column = (seq_length - seq_ch_taps + 1 );  s_row    =  seq_max_ant * seq_ch_taps;  par    = training_seq[ type ].shsish;  seq    = training_seq[ type ].sequence;  for (i = 0; i < s_row; i++) {    double_real = 0.;    double_imag = 0.;    for (j = 0; j < s_column; j++) {      double_real += par[ 2*j + 2*i*s_column ] * mid[ -j ].real - 	par[ 2*j + 1 + 2*i*s_column ] * mid[ -j ].imag;      double_imag += par[ 2*j + 2*i*s_column ] * mid[ -j ].imag + 	par[ 2*j + 1 + 2*i*s_column ] * mid[ -j ].real;    }    channel[ i ].real = double_real; //Quantization double -> int    channel[ i ].imag = double_imag;  }  // Compute the matched filter: zero-forcing equalisation  // channel_inv = ifft( 1./ fft( channel ) )  // We have to do this once for every received channel  {    f_len = pow( 2, floor( logb( seq_ch_taps ) ) );    for ( ch=0; ch<num_of_antennas; ch++ ){      int ch_off = ch * seq_ch_taps;      complex double *f_inv = f_invs[ch];      if ( out[ ch ] ) {	f_abs = 0;	for ( i=0; i<f_len; i++ ){	  f_inv[i] = SC_TO_CD( channel[ i+ch_off ] );