chest_rcv_mimo.c

软件无线电的平台

/***************************************************************************
 *   chest_rcv.c  - MIMO channel estimation module 
 *   begin          :  03/08/05
 *   Author         : Selvavinayagam Gunabalan
 *   emails         : [e-mail]
 ***************************************************************************/
/***************************************************************************
 *                                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 4


typedef 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 ] );