decode.c

剑桥大学David J.C. MacKay 个人网站公布的2006年的代码

/* DECODE.C - Decode blocks of received data. */

/* Copyright (c) 2000, 2001 by Radford M. Neal 
 *
 * Permission is granted for anyone to copy, use, modify, or distribute this
 * program and accompanying programs and documents for any purpose, provided 
 * this copyright notice is retained and prominently displayed, along with
 * a note saying that the original programs are available from Radford Neal's
 * web page, and note is made of any changes made to the programs.  The
 * programs and documents are distributed without any warranty, express or
 * implied.  As the programs were written for research purposes only, they have
 * not been tested to the degree that would be advisable in any important
 * application.  All use of these programs is entirely at the user's own risk.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "rand.h"
#include "alloc.h"
#include "blockio.h"
#include "open.h"
#include "mod2sparse.h"
#include "mod2dense.h"
#include "mod2convert.h"
#include "channel.h"
#include "rcode.h"
#include "check.h"
#include "dec.h"


void usage(void);


/* MAIN PROGRAM. */

int main
( int argc,
  char **argv
)
{
  char *pchk_file, *rfile, *dfile, *pfile;
  char **meth;
  FILE *rf, *df, *pf;

  char *dblk, *pchk;
  double *lratio;
  double *bitpr;

  double *awn_data;		/* Places to store channel data */
  int *bsc_data;

  unsigned iters;		/* Unsigned because can be huge for enum */
  double tot_iter;		/* Double because can be huge for enum */
  double chngd, tot_changed;	/* Double because can be fraction if lratio==1*/

  int tot_valid;
  char junk;
  int valid;

  int i, j, k;

  /* Look at arguments up to the decoding method specification. */

  table = 0;
  if (argc>1 && strcmp(argv[1],"-t")==0)
  { table = 1;
    argc -= 1;
    argv += 1;
  }
  if (argc>1 && strcmp(argv[1],"-T")==0)
  { table = 2;
    argc -= 1;
    argv += 1;
  }

  if (!(pchk_file = argv[1])
   || !(rfile = argv[2])
   || !(dfile = argv[3]))
  { usage();
  }

  if (argv[4]==0 || argv[5]==0) usage();

  k = channel_parse(argv+4,argc-4);
  if (k<=0)
  { pfile = argv[4];
    k = channel_parse(argv+5,argc-5);
    if (k<=0) usage();
    meth = argv+5+k;
  }
  else
  { pfile = 0;
    meth = argv+4+k;
  }

  /* Look at the specification of the decoding method, which starts at meth and
     continues to the end of the command line (marked by a zero pointer). */

  if (!meth[0]) usage();

  if (strcmp(meth[0],"prprp")==0)
  { dec_method = Prprp;
    if (!meth[1] || sscanf(meth[1],"%d%c",&max_iter,&junk)!=1 || meth[2]) 
    { usage();
    }
  }
  else if (strcmp(meth[0],"enum-block")==0)
  { dec_method = Enum_block;
    if (!(gen_file = meth[1]) || meth[2]) usage();
  }
  else if (strcmp(meth[0],"enum-bit")==0)
  { dec_method = Enum_bit;
    if (!(gen_file = meth[1]) || meth[2]) usage();
  }
  else 
  { usage();
  }

  /* Check that we aren't overusing standard input or output. */

  if ((strcmp(pchk_file,"-")==0) 
    + (strcmp(rfile,"-")==0) > 1)
  { fprintf(stderr,"Can't read more than one stream from standard input\n");
    exit(1);
  }
  if ((table>0) 
    + (strcmp(dfile,"-")==0) 
    + (pfile!=0 && strcmp(pfile,"-")==0) > 1)
  { fprintf(stderr,"Can't send more than one stream to standard output\n");
    exit(1);
  }

  /* Read parity check file. */

  read_pchk(pchk_file);

  if (N<=M)
  { fprintf(stderr,
     "Number of bits (%d) should be greater than number of checks (%d)\n",N,M);
    exit(1);
  }

  /* Open file of received data. */

  rf = open_file_std(rfile,"r");
  if (rf==NULL)
  { fprintf(stderr,"Can't open file of received data: %s\n",rfile);
    exit(1);
  }

  /* Create file for decoded data. */

  df = open_file_std(dfile,"w");
  if (df==NULL)
  { fprintf(stderr,"Can't create file for decoded data: %s\n",dfile);
    exit(1);
  }

  /* Create file for bit probabilities, if specified. */

  if (pfile)
  { pf = open_file_std(pfile,"w");
    if (pf==NULL)
    { fprintf(stderr,"Can't create file for bit probabilities: %s\n",pfile);
      exit(1);
    }
  }

  /* Allocate space for data from channel. */

  switch (channel)
  { case BSC:
    { bsc_data = chk_alloc (N, sizeof *bsc_data);
      break;
    }
    case AWGN: case AWLN:
    { awn_data = chk_alloc (N, sizeof *awn_data);
      break;
    }
    default:
    { abort();
    }
  }

  /* Allocate other space. */

  dblk   = chk_alloc (N, sizeof *dblk);
  lratio = chk_alloc (N, sizeof *lratio);
  pchk   = chk_alloc (M, sizeof *pchk);
  bitpr  = chk_alloc (N, sizeof *bitpr);

  /* Print header for summary table. */

  if (table==1)
  { printf("  block iterations valid  changed\n");
  }

  /* Do the setup for the decoding method. */

  switch (dec_method)
  { case Prprp: 
    { prprp_decode_setup();
      break;
    }
    case Enum_block: case Enum_bit:
    { enum_decode_setup();
      break;
    }
    default: abort();
  }

  /* Read received blocks, decode, and write decoded blocks. */

  tot_iter = 0;
  tot_valid = 0;
  tot_changed = 0;

  for (block_no = 0; ; block_no++)
  { 
    /* Read block from received file, exit if end-of-file encountered. */

    for (i = 0; i<N; i++)
    { int c;
      switch (channel)
      { case BSC:  
        { c = fscanf(rf,"%1d",&bsc_data[i]); 
          break;
        }
        case AWGN: case AWLN:
        { c = fscanf(rf,"%lf",&awn_data[i]); 
          break;
        }
        default: abort();
      }
      if (c==EOF) 
      { if (i>0)
        { fprintf(stderr,
          "Warning: Short block (%d long) at end of received file ignored\n",i);
        }
        goto done;
      }
      if (c<1 || channel==BSC && bsc_data[i]!=0 && bsc_data[i]!=1)
      { fprintf(stderr,"File of received data is garbled\n");
        exit(1);
      }
    }

    /* Find likelihood ratio for each bit. */

    switch (channel)
    { case BSC:
      { for (i = 0; i<N; i++)
        { lratio[i] = bsc_data[i]==1 ? (1-error_prob) / error_prob
                                     : error_prob / (1-error_prob);
        }
        break;
      }
      case AWGN:
      { for (i = 0; i<N; i++)
        { lratio[i] = exp(2*awn_data[i]/(std_dev*std_dev));
        }
        break;
      }
      case AWLN:
      { for (i = 0; i<N; i++)
        { double e, d1, d0;
          e = exp(-(awn_data[i]-1)/lwidth);
          d1 = 1 / ((1+e)*(1+1/e));
          e = exp(-(awn_data[i]+1)/lwidth);
          d0 = 1 / ((1+e)*(1+1/e));
          lratio[i] = d1/d0;
        }
        break;
      }
      default: abort();
    }

    /* Try to decode using the specified method. */

    switch (dec_method)
    { case Prprp:
      { iters = prprp_decode (H, lratio, dblk, pchk, bitpr);
        break;
      }
      case Enum_block: case Enum_bit:
      { iters = enum_decode (lratio, dblk, bitpr, dec_method==Enum_block);
        break;
      }
      default: abort();
    }

    /* See if it worked, and how many bits were changed. */

    valid = check(H,dblk,pchk)==0;

    chngd = changed(lratio,dblk,N);

    tot_iter += iters;
    tot_valid += valid;
    tot_changed += chngd;

    /* Print summary table entry. */

    if (table==1)
    { printf ("%7d %10f    %d  %8.1f\n",
        block_no, (double)iters, valid, (double)chngd);
        /* iters is printed as a double to avoid problems if it's >= 2^31 */
      fflush(stdout);
    }

    /* Write decoded block. */

    blockio_write(df,dblk,N);

    /* Write bit probabilities, if asked to. */

    if (pfile)
    { for (j = 0; j<N; j++)
      { fprintf(pf," %.5f",bitpr[j]);
      }
      fprintf(pf,"\n");
    }
  }

  /* Finish up. */

done: 
  fprintf(stderr,
  "Decoded %d blocks, %d valid.  Average %.1f iterations, %.0f%% bit changes\n",
   block_no, tot_valid, (double)tot_iter/block_no, 
   100.0*(double)tot_changed/(N*block_no));

  if (ferror(df) || fclose(df)!=0)
  { fprintf(stderr,"Error writing decoded blocks to %s\n",dfile);
    exit(1);
  }

  if (pfile)
  { if (ferror(pf) || fclose(pf)!=0)
    { fprintf(stderr,"Error writing bit probabilities to %s\n",dfile);
      exit(1);
    }
  }

  exit(0);
}


/* PRINT USAGE MESSAGE AND EXIT. */


void usage(void)
{ fprintf(stderr,"Usage:\n");
  fprintf(stderr,
"  decode [ -t ] pchk-file received-file decoded-file [ bp-file ] channel method\n");
  channel_usage();
  fprintf(stderr,
"Method:  enum-block gen-file | enum-bit gen-file | prprp [-]max-iterations\n");
  exit(1);
}