ibn.c

快速傅立叶变换程序代码,学信号的同学,可要注意了

From utstat.toronto.edu!radford Sun Oct  8 19:41:53 1995
Subject: Sucessful decoding above R0
To: [e-mail]
From: "Radford Neal" <[e-mail]>
Cc: [e-mail] (Radford Neal)

David,

I've generalized my version of the infinite belief net simulator, so
that I can do try decoding when f_s is not equal to f_n (as well as
various other tricks).

I've just played a little so far, but have some interesting results
(assuming there are no bugs, of course).  I thought I'd try the
approach of fixing the noise level and the desired information rate,
and seeing whether the code can be fiddled to achieve error-free
transmission.  I've apparently succeeded with f_n = 0.01 and an
information rate of 0.79.  Note that the capacity with f_n = 0.01 is
0.9192 and R0 is 0.7382.

The code that does this is just the t=3, tr = 6 one, but with
f_s=0.2370.  I assume that the 6 bits that go into each check are
evenly divided into 3 noise bits and 3 signal bits.  This seems like
it could be arranged (near enough).

I've appended my new program, and the results of the run.  The results
are for a sample size of 20,000; a sample size of 10,000 gives similar
results.  Let me know if you see any bugs!  

   Radford

-------------------------------------------------------------------------------

/* IBN.C - Simulate probability propagation in infinite belief network. */

/* Usage: 

   ibn [ -p ] seed noise-pr [ sig-pr/rate ] noise-bits [ sig-bits [ sig-chks ] ]

   Simulates belief-net decoding with an infinite network.  The output of the 
   program shows the progress of decoding at each iteration, up to I iterations
   (see below).  If the optional "-p" argument is included, the output is shorn 
   of headings, so as to be suitable for plotting, and gives only the iteration
   number and the probability of an error in decoding a signal bit after that 
   number of iterations.

   The computations are done using a Monte Carlo technique with sample size
   S (see below), using the given random number seed.  

   The noise probability is given by the following argument, and may be
   followed by another argument giving the probability of a 1 in the signal, 
   or if it begins with a minus sign, the desired information transmission 
   rate (symbol rate times entropy).  If this argument is omitted, the
   default is for the signal to have the same properties as the noise.

   The next argument gives the number of noise bits participating in each
   check.  A value of one gives a simple sparse matrix code; values above one 
   give standard MN codes.  The number of checks applying to each noise bit 
   is assumed to be the same as this.

   The number of signal bits participating in each check follows, with the
   default being the same as the number of noise bits.  The number of 
   checks that each signal bit is part of follows, with the default again being
   the same as the number of noise bits in each check.  The symbol rate of 
   the code is sig-bits/sig-chks.
*/

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

#include "rand.h"


/* ADJUSTABLE PARAMETERS. */

#define I 35		/* Number of iterations */
#define S 20000		/* Sample size */


/* STATE OF THE SIMULATION. */

/* The decoder's probability for whether a bit is a 1 given particular checks 
   farther out in the tree, for a sample of actual checks computed from 
   this bit and bits farther out.  Four samples of size S are maintained:
   for signal and noise bits, when the bit in question is actually a 0 or 
   actually a 1. */

float np[2][S], sp[2][S];

/* The decoder's likelihood ratio for a check having its actual value given
   that its inner parent bit is a 1 vs. a 0, for a sample of actual checks
   computed from this bit and bits farther out.  Four samples of size S are 
   maintained: for when the inner bit is a signal or noise bit, when this bit 
   is is actually a 0 or actually a 1. */

float nl[2][S], sl[2][S];


/* ENTROPY FUNCTION (IN BITS). */

#define H(p) \
  (((p)==0 || 1-(p)==0 ? 0 : -(p)*log(p)-(1-(p))*log(1-(p))) / log(2.0))


/* PRINT USAGE MESSAGE AND EXIT. */

static void usage()
{
  fprintf(stderr,"Usage:\n");
  fprintf(stderr,
   "  ibn [ -p ] seed noise-pr [ sig-pr/rate ] noise-bits [ sig-bits [ sig-chks ] ]\n");
  exit(1);
}


/* MAIN PROGRAM. */

void main
( int argc,
  char **argv
)
{
  double noise_pr, sig_pr, sym_rate, info_rate;
  int noise_chks, noise_bits, sig_chks, sig_bits;

  int plot, seed;

  double h[2], e[2];
  double l, p, p1;
  int chk, bit;

  int i, b, s, n, m;

  /* Get arguments from command line. */

  i = 1;

  plot = 0;
  if (i<argc && strcmp(argv[i],"-p")==0) 
  { plot = 1;
    i += 1;
  }

  if (i>=argc || (seed = atoi(argv[i++]))<=0) usage();

  if (i>=argc || (noise_pr = atof(argv[i++]))<=0 || noise_pr>=1) usage();

  sig_pr = noise_pr;
  if (i<argc && strchr(argv[i],'.'))
  { sig_pr = atof(argv[i++]);
    if (sig_pr==0 || sig_pr>=1 || sig_pr<=-1) usage();
  }

  if (i>=argc || (noise_bits = atoi(argv[i++]))<=0) usage();

  noise_chks = noise_bits;

  sig_bits = noise_bits;
  sig_chks = noise_bits;

  if (i<argc)
  { if ((sig_bits = atoi(argv[i++]))<=0) usage();
  }

  if (i<argc)
  { if ((sig_chks = atoi(argv[i++]))<=0) usage();
  }

  if (i!=argc) usage();

  /* Figure out symbol rate and find or figure out information rate. */

  sym_rate = (double) sig_bits / sig_chks;

  if (sig_pr<0)
  { double nh, lp, hp, mp;
    nh = (-sig_pr) / sym_rate;
    if (nh>1)
    { fprintf(stderr,
       "Requested information rate of %.4lf exceeds the symbol rate of %.4lf\n",
        -sig_pr, sym_rate);
      exit(1);
    }
    lp = 0; hp = 0.5;
    while (hp-lp>1e-10)
    { mp = (hp+lp)/2;
      if (H(mp)<nh) lp = mp;
      else hp = mp;
    }
    sig_pr = (lp+hp)/2;
  }

  info_rate = sym_rate * H(sig_pr);
  
  /* Print information, unless producing plot output. */

  if (!plot)
  { printf(
    "\nNoise-pr = %.4lf  Sig-pr = %.4lf  Sym-rate = %.4lf  Info-rate = %.4lf\n",
     noise_pr, sig_pr, sym_rate, info_rate);
    printf(
    "\nNoise-bits = Noise-chks = %d  Sig-bits = %d  Sig-chks = %d\n",
     noise_bits, sig_bits, sig_chks);
    printf(
    "\nCapacity = %.4lf  R0 = %.4lf  GV = %.4lf\n",
     1-H(noise_pr), 1-log(1+2*sqrt(noise_pr*(1-noise_pr)))/log(2.0),
     noise_pr<0.25 ? 1-H(2*noise_pr) : 0);
    printf(
    "\nSeed = %d  Sample-size = %d  Iterations = %d\n\n", 
     seed, S, I);
  }

  /* Set up probabilities for outer fringe of bits - the decoder just
     knows their prior probability of being 1. */

  for (b = 0; b<2; b++)
  { for (s = 0; s<S; s++) 
    { np[b][s] = noise_pr;
      sp[b][s] = sig_pr;
    }
  }

  /* Propagate inwards toward signal/noise bit at root of tree. */

  rand_seed(seed);

  for (i = 1; i<=I; i++)
  { 
    /* Find samples for likelihood ratios of randomly generated checks. */

    for (b = 0; b<2; b++)
    { for (s = 0; s<S; s++)
      { 
        /* When inner bit is a noise bit. */

        chk = b;
        p1 = 0;
        for (n = 0; n<noise_bits-1; n++)
        { bit = rand_uniform()<noise_pr;
          chk ^= bit;
          p = np[bit][rand_int(S)];
          p1 = p1*(1-p) + (1-p1)*p;
        }
        for (n = 0; n<sig_bits; n++)
        { bit = rand_uniform()<sig_pr;
          chk ^= bit;
          p = sp[bit][rand_int(S)];
          p1 = p1*(1-p) + (1-p1)*p;
        }
        nl[b][s] = chk ? (1-p1)/p1 : p1/(1-p1);

        /* When inner bit is a signal bit. */

        chk = b;
        p1 = 0;
        for (n = 0; n<noise_bits; n++)
        { bit = rand_uniform()<noise_pr;
          chk ^= bit;
          p = np[bit][rand_int(S)];
          p1 = p1*(1-p) + (1-p1)*p;
        }
        for (n = 0; n<sig_bits-1; n++)
        { bit = rand_uniform()<sig_pr;
          chk ^= bit;
          p = sp[bit][rand_int(S)];
          p1 = p1*(1-p) + (1-p1)*p;
        }
        sl[b][s] = chk ? (1-p1)/p1 : p1/(1-p1);
      }
    }

    /* Find samples for decoder's probabilities for bits. */

    for (b = 0; b<2; b++)
    { for (s = 0; s<S; s++)
      { 
        /* For noise bits. */

        l = noise_pr/(1-noise_pr);
        for (m = 0; m<noise_chks-1; m++)
        { l *= nl[b][rand_int(S)];
        }
        np[b][s] = finite(l) ? l/(1+l) : 1;

        /* For signal bits. */

        l = sig_pr/(1-sig_pr);
        for (m = 0; m<sig_chks-1; m++)
        { l *= sl[b][rand_int(S)];
        }
        sp[b][s] = finite(l) ? l/(1+l) : 1;
      }
    }

    /* Compute average entropy and probability of decoding errors for signal
       bits, separately for bits that are actually 0 or 1, and averaged. */

    for (b = 0; b<2; b++)
    { h[b] = 0; e[b] = 0;
      for (s = 0; s<S; s++)  
      { h[b] += H(sp[b][s]);
        e[b] += (sp[b][s]>=0.5) != b;
      }
      h[b] /= S; e[b] /= S;
    }

    /* Print results for this iteration.  When plotting, print only averaged 
       probability of error. */

    if (plot)
    { printf ("%4d %.5lf\n", i, sig_pr*e[1] + (1-sig_pr)*e[0]);
    }
    else
    { printf 
       ("%4d - h0: %.5lf h1: %.5lf h: %.5lf - e0: %.5lf e1: %.5lf e: %.5lf\n",
         i, h[0], h[1], sig_pr*h[1] + (1-sig_pr)*h[0],
         e[0], e[1], sig_pr*e[1] + (1-sig_pr)*e[0]);
    }
  }

  exit(0);
}

-------------------------------------------------------------------------------

Output of the command "ibn 1 0.01 -0.79 3":

Noise-pr = 0.0100  Sig-pr = 0.2370  Sym-rate = 1.0000  Info-rate = 0.7900

Noise-bits = Noise-chks = 3  Sig-bits = 3  Sig-chks = 3

Capacity = 0.9192  R0 = 0.7382  GV = 0.8586

Seed = 1  Sample-size = 20000  Iterations = 35

   1 - h0: 0.68641 h1: 0.82742 h: 0.71983 - e0: 0.00000 e1: 1.00000 e: 0.23699
   2 - h0: 0.63510 h1: 0.82160 h: 0.67930 - e0: 0.04830 e1: 0.75330 e: 0.21538
   3 - h0: 0.58933 h1: 0.81403 h: 0.64259 - e0: 0.05175 e1: 0.70080 e: 0.20557
   4 - h0: 0.56741 h1: 0.80427 h: 0.62354 - e0: 0.06130 e1: 0.62630 e: 0.19520
   5 - h0: 0.53971 h1: 0.79577 h: 0.60039 - e0: 0.06090 e1: 0.60520 e: 0.18989
   6 - h0: 0.52358 h1: 0.78431 h: 0.58537 - e0: 0.06425 e1: 0.56740 e: 0.18349
   7 - h0: 0.50216 h1: 0.77253 h: 0.56623 - e0: 0.06660 e1: 0.55220 e: 0.18168
   8 - h0: 0.49639 h1: 0.76103 h: 0.55911 - e0: 0.06745 e1: 0.51880 e: 0.17442
   9 - h0: 0.48216 h1: 0.75617 h: 0.54710 - e0: 0.06510 e1: 0.50540 e: 0.16945
  10 - h0: 0.47093 h1: 0.74460 h: 0.53579 - e0: 0.06675 e1: 0.48995 e: 0.16704
  11 - h0: 0.45511 h1: 0.73873 h: 0.52233 - e0: 0.06390 e1: 0.46565 e: 0.15911
  12 - h0: 0.43997 h1: 0.72688 h: 0.50797 - e0: 0.05790 e1: 0.44505 e: 0.14965
  13 - h0: 0.41255 h1: 0.71595 h: 0.48445 - e0: 0.05880 e1: 0.42345 e: 0.14522
  14 - h0: 0.40265 h1: 0.69344 h: 0.47156 - e0: 0.06550 e1: 0.40005 e: 0.14479
  15 - h0: 0.38945 h1: 0.67890 h: 0.45804 - e0: 0.06345 e1: 0.38495 e: 0.13964
  16 - h0: 0.37187 h1: 0.66575 h: 0.44152 - e0: 0.05515 e1: 0.37755 e: 0.13156
  17 - h0: 0.35886 h1: 0.64845 h: 0.42749 - e0: 0.05740 e1: 0.34920 e: 0.12655
  18 - h0: 0.33859 h1: 0.63008 h: 0.40767 - e0: 0.05640 e1: 0.32955 e: 0.12113
  19 - h0: 0.32032 h1: 0.61046 h: 0.38908 - e0: 0.05480 e1: 0.30415 e: 0.11389
  20 - h0: 0.30517 h1: 0.58263 h: 0.37093 - e0: 0.04975 e1: 0.28035 e: 0.10440
  21 - h0: 0.28010 h1: 0.55664 h: 0.34564 - e0: 0.04645 e1: 0.25700 e: 0.09635
  22 - h0: 0.25407 h1: 0.51551 h: 0.31603 - e0: 0.04760 e1: 0.22875 e: 0.09053
  23 - h0: 0.23150 h1: 0.47277 h: 0.28868 - e0: 0.04445 e1: 0.20355 e: 0.08216
  24 - h0: 0.20224 h1: 0.43220 h: 0.25674 - e0: 0.03905 e1: 0.17985 e: 0.07242
  25 - h0: 0.17741 h1: 0.37209 h: 0.22354 - e0: 0.03655 e1: 0.15020 e: 0.06348
  26 - h0: 0.14176 h1: 0.32035 h: 0.18408 - e0: 0.02875 e1: 0.12350 e: 0.05120
  27 - h0: 0.10537 h1: 0.25183 h: 0.14008 - e0: 0.02090 e1: 0.09175 e: 0.03769
  28 - h0: 0.06891 h1: 0.17422 h: 0.09387 - e0: 0.01365 e1: 0.05700 e: 0.02392
  29 - h0: 0.03257 h1: 0.09555 h: 0.04750 - e0: 0.00720 e1: 0.03035 e: 0.01269
  30 - h0: 0.01226 h1: 0.03152 h: 0.01682 - e0: 0.00285 e1: 0.00920 e: 0.00435
  31 - h0: 0.00178 h1: 0.00546 h: 0.00265 - e0: 0.00025 e1: 0.00135 e: 0.00051
  32 - h0: 0.00016 h1: 0.00021 h: 0.00017 - e0: 0.00000 e1: 0.00000 e: 0.00000
  33 - h0: 0.00000 h1: 0.00000 h: 0.00000 - e0: 0.00000 e1: 0.00000 e: 0.00000
  34 - h0: 0.00000 h1: 0.00000 h: 0.00000 - e0: 0.00000 e1: 0.00000 e: 0.00000
  35 - h0: 0.00000 h1: 0.00000 h: 0.00000 - e0: 0.00000 e1: 0.00000 e: 0.00000