加载中…
个人资料
我孙子的爷爷
我孙子的爷爷
  • 博客等级:
  • 博客积分:0
  • 博客访问:27,747
  • 关注人气:1
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
正文 字体大小:

仿真卷积码和Viterbi解码的C程序

(2006-08-10 19:50:28)

// ------------------------------------------------------------------------
// File: viterbi_binary_nsc.c
// Date: April 1, 2002
// Description: Viterbi decoder, maximum likelihood decoding. For a
//              rate-1/n convolutional code. Correlation metric.
//
// NOTE: The generator polynomials are given in hexadecimal and in
//       reversed order. For example, the rate-1/2 standard code
//       has generator polynomials (octal) (g1,g2)=(133,171). Then
//
//           Base 8       Reverse      Base 16
//       133 = 1011011     --->     1101101 = 6d
//       171 = 1111001     --->     1001111 = 4f
//
//       These values can be found in file "trellis_64states.data"
//
//       Also, the order in the file is (g2,g1). That is, the last
//       generator in the first row, down to the first generator in
//       the last row.
//
//       With this format, the software implementation of the
//       algorithm is simplified.
//
// NOTE: Maximum of 1024 states; max. truncation length = 1024
//
// ------------------------------------------------------------------------
 






#include "stdio.h"
#include "math.h"
#include "float.h"
#include "limits.h"
#include "time.h" 
 
#define MAX_RANDOM LONG_MAX
 
// Use the compiling directive -DSHOW_PROGRESS to see how the decoder
// converges to the decoded sequence by monitoring the survivor memory
#ifdef SHOW_PROGRESS
#define DELTA 1
#endif
 
int k2=1, n2, m2;
int NUM_STATES, OUT_SYM, NUM_TRANS;
long TRUNC_LENGTH;
 
double RATE;
double INIT_SNR, FINAL_SNR, SNR_INC;
long NUMSIM;
 
FILE *fp, *fp_ber;
 
int g2[10][10];
unsigned int memory2, output;           
unsigned int data2;                     
 
unsigned long seed;                      
unsigned int data_symbol[1024];         
long index;                             
double transmitted[10];                 
double snr, amp;
double received[10];                    
 
int fflush();
 
// Data structures used for trellis sections and survivors
struct trel {
  int init;               
  int data;               
  int final;              
  double output[10]; 
};
struct surv {
  double metric;          
  int data[1024]; 
  int state[1024]; 
};
 
// A trellis section is an array of branches, indexed by an initial
//   state and a k_2-bit input data. The values read
//   are the final state and the output symbols
struct trel trellis[1024][100];
 

struct surv survivor[1024], surv_temp[1024];
 

void encoder2(void);          
int random_data(void);        
void transmit(void);          
void awgn(void);              
void viterbi(void);           
double comp_metric(double rec, double ref);
void open_files(void);        
void close_files(void);       
 
 
main(int argc, char *argv[])
{
 
register int i, j, k;
int init, data, final, output;
register int error;
unsigned long error_count;
char name1[40], name2[40];
 
  // Command line processing
  if (argc != 8)
    {
      printf("Usage %s file_input file_output truncation snr_init snr_final snr_inc num_sim\n", argv[0]);
      exit(0);
    }
 
  sscanf(argv[1],"%s", name1);
  sscanf(argv[2],"%s", name2);
  sscanf(argv[3],"%ld", &TRUNC_LENGTH);
  sscanf(argv[4],"%lf", &INIT_SNR);
  sscanf(argv[5],"%lf", &FINAL_SNR);
  sscanf(argv[6],"%lf", &SNR_INC);
  sscanf(argv[7],"%ld", &NUMSIM);
 
printf("\nSimulation of Viterbi decoding over an AWGN channel\n");
printf("%ld simulations per Eb/No (dB) point\n", NUMSIM);
 
 
  fp_ber = fopen(name2,"w");
 
 
  if (!(fp = fopen(name1,"r")))
    {
    printf("Error opening file!\n");
    exit(1);
    }
 
  fscanf(fp, "%d %d", &n2, &m2);
 
  RATE = 1.0 / (double) n2;
 
  fscanf(fp, "%d %d %d", &NUM_STATES, &OUT_SYM, &NUM_TRANS);
  for (j=0; j
    fscanf(fp, "%x", &g2[j][0]);
 
  printf("\n%d-state rate-1/%d binary convolutional encoder\n",
          NUM_STATES, n2);
  printf("with generator polynomials ");
  for (j=0; j
  printf("\n\nDecoding depth = %ld\n\n", TRUNC_LENGTH);
 
 
 
  for (j=0; j
    for (k=0; k
      {
       
        fscanf(fp,"%d %d %d", &trellis[j][k].init, &trellis[j][k].data,
        &trellis[j][k].final);
       
        for (i=0; i < OUT_SYM; i++)
          {
             fscanf(fp,"%d", &data);
             trellis[j][k].output[i] = (double) data;
          }
      }
     
 
  fclose(fp);
 
#ifdef PRINT
  for (j=0; j  
    for (k=0; k 
      {
        printf("= = =  ", trellis[j][k].init,
        trellis[j][k].data, trellis[j][k].final);
        for (i=0; i < OUT_SYM; i++)
          printf("%2.0lf ", trellis[j][k].output[i]);
        printf("\n");
       }
     
#endif
 
 
 
  snr = INIT_SNR;
 
 
 
  while ( snr < (FINAL_SNR+1.0e-6) )
    {
      amp = sqrt( 2.0 * RATE * pow(10.0, (snr/10.0)) );
 
     
      time(&seed);
      srandom(seed);
     
     
     
      for (i=0; i
        data_symbol[i]=0;
     
     
     
      for (i=0; i
        {
          survivor[i].metric = 0.0;            
          for (j=0; j
            {
              survivor[i].data[j] = 0;       
              survivor[i].state[j] = 0;      
            }
        }
     
 
     
      index = 0;
     
     
      memory2 = 0;
 
      
      error_count = 0;
 
 
     
     
      while (index < NUMSIM)
        {
       
         
          data_symbol[index % TRUNC_LENGTH] = random_data();
         
         
          transmit();
         
         
          awgn();
 
         
          viterbi();
          
          index += 1;          
         
         
          error = survivor[0].data[index % TRUNC_LENGTH]
                                      ^ data_symbol[index % TRUNC_LENGTH];
          if (error)
            error_count++;
         
#ifdef SHOW_PROGRESS
          if ( (index % DELTA) == 0 )
            {
              for (j=0; j
                {
                  printf("= -   ", (index % TRUNC_LENGTH), j);
                  for (i=0; i
                    printf("%x", survivor[j].data[i]);
                  printf(" %ld\n", error_count);
                }
            }
#endif
        }
     
      printf("%f  .4e\n",snr, ((double) error_count/(double) index) );
     
      fprintf(fp_ber, "%f .4e\n", snr, ((double)error_count /(double)index));
             
 
      fflush(stdout);
      fflush(fp_ber);
 
      snr += SNR_INC;
     
    }
 
  fclose(fp_ber);
}
 
 
 
 
 
int random_data()
{
 
  return( (random() >> 5) & 1 );
}
 
 
 
 
 
 
 
void transmit()
{
 
  int i;
 
  encoder2();
 
 
  for (i=(OUT_SYM-1); i >=0; i--)
    if ( (output >> i) & 0x01 )
      transmitted[OUT_SYM-1-i] = -1.0; 
    else
      transmitted[OUT_SYM-1-i] =  1.0; 
}
 
 
 
 
 
 
void encoder2()
{
  // Binary convolutional encoder, rate 1/n2
  register int i, j, result, temp;
 
  temp = memory2;
  output = 0;
 
  temp = (temp<<1) ^ data_symbol[index % TRUNC_LENGTH];
 
  for (i=0; i
   {
     result = 0;
     for (j=m2; j>=0; j--)
       result ^= ( ( temp & g2[i][0] ) >> j ) & 1;
     output = ( output<<1 ) ^ result;
   }
 
  memory2 = temp ;
 
}
 
 
 
 
 
 
void awgn()
{
 
  double u1,u2,s,noise,randmum;
  int i;
 
  for (i=0; i
    {
 
      do
        {      
          randmum = (double)(random())/MAX_RANDOM;
          u1 = randmum*2.0 - 1.0;
          randmum = (double)(random())/MAX_RANDOM;
          u2 = randmum*2.0 - 1.0;
          s = u1*u1 + u2*u2;
        } while( s >= 1);
      noise = u1 * sqrt( (-2.0*log(s))/s );
      received[i] = transmitted[i] + noise/amp;
    }
}
 
 
 
 
 
 
 
 
void viterbi()
{
  double aux_metric, surv_metric[NUM_STATES];
  register int i,j,k;
 
  for (i=0; i
    surv_metric[i] = -DBL_MAX;
 
  for (i=0; i
    {
      for (j=0; j
        {
 
         
          aux_metric = 0.0;
          for (k=(OUT_SYM-1); k>=0; k--)
            aux_metric += comp_metric(received[k],trellis[i][j].output[k]);
          aux_metric += survivor[i].metric;
 
         
         
 
          if ( aux_metric > surv_metric[trellis[i][j].final] )
            {
 
              surv_metric[trellis[i][j].final] = aux_metric;
 
             
              for (k=0; k
                {
                  surv_temp[trellis[i][j].final].data[k] =
                    survivor[i].data[k];
                  surv_temp[trellis[i][j].final].state[k] =
                    survivor[i].state[k];
                }
                  surv_temp[trellis[i][j].final].data[index%TRUNC_LENGTH] = j;
                  surv_temp[trellis[i][j].final].state[index%TRUNC_LENGTH] =
                    trellis[i][j].final;
 
            }
        }
    }
 
  for (i=0; i
    {
      
      survivor[i].metric = surv_metric[i];
      for (k=0; k
        {
          survivor[i].data[k] = surv_temp[i].data[k];
          survivor[i].state[k] = surv_temp[i].state[k];
        }
    }
 
}
 
 
 
 
 
 
 
double comp_metric(double rec, double ref)
{
 


 
  if (ref == 1.0)
    return(rec);
  else
    return(-rec);
}
 
 

0

阅读 收藏 喜欢 打印举报/Report
  

新浪BLOG意见反馈留言板 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

新浪公司 版权所有