Imported Upstream version 3.2.2
[debian/gnuradio] / gnuradio-core / src / lib / viterbi / metrics.c
1 /*
2  * Copyright 1995 Phil Karn, KA9Q
3  * Copyright 2008 Free Software Foundation, Inc.
4  * 
5  * This file is part of GNU Radio
6  * 
7  * GNU Radio is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 3, or (at your option)
10  * any later version.
11  * 
12  * GNU Radio is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with GNU Radio; see the file COPYING.  If not, write to
19  * the Free Software Foundation, Inc., 51 Franklin Street,
20  * Boston, MA 02110-1301, USA.
21  */
22
23 /* 
24  * Generate metric tables for a soft-decision convolutional decoder
25  * assuming gaussian noise on a PSK channel.
26  *
27  * Works from "first principles" by evaluating the normal probability
28  * function and then computing the log-likelihood function
29  * for every possible received symbol value
30  *
31  */
32
33 #ifdef HAVE_CONFIG_H
34 #include <config.h>
35 #endif
36
37 /* Symbols are offset-binary, with 128 corresponding to an erased (no
38  * information) symbol
39  */
40 #define OFFSET  128
41
42 #include <stdlib.h>
43 #include <math.h>
44
45 /* Normal function integrated from -Inf to x. Range: 0-1 */
46 #define normal(x)       (0.5 + 0.5*erf((x)/M_SQRT2))
47
48 /* Logarithm base 2 */ 
49 #define gr_log2(x) (log(x)*M_LOG2E)
50
51 /* Generate log-likelihood metrics for 8-bit soft quantized channel
52  * assuming AWGN and BPSK
53  */
54 void
55 gen_met(int mettab[2][256],     /* Metric table, [sent sym][rx symbol] */
56         int amp,                /* Signal amplitude, units */
57         double esn0,            /* Es/N0 ratio in dB */
58         double bias,            /* Metric bias; 0 for viterbi, rate for sequential */
59         int scale)              /* Scale factor */
60 {
61   double noise;
62   int s,bit;
63   double metrics[2][256];
64   double p0,p1;
65   
66   /* Es/N0 as power ratio */
67   esn0 = pow(10.,esn0/10);
68   
69   noise = 0.5/esn0;     /* only half the noise for BPSK */
70   noise = sqrt(noise);  /* noise/signal Voltage ratio */
71   
72   /* Zero is a special value, since this sample includes all
73    * lower samples that were clipped to this value, i.e., it
74    * takes the whole lower tail of the curve 
75    */
76   p1 = normal(((0-OFFSET+0.5)/amp - 1)/noise);  /* P(s|1) */
77   
78   /* Prob of this value occurring for a 0-bit */        /* P(s|0) */
79   p0 = normal(((0-OFFSET+0.5)/amp + 1)/noise);
80   metrics[0][0] = gr_log2(2*p0/(p1+p0)) - bias;
81   metrics[1][0] = gr_log2(2*p1/(p1+p0)) - bias;
82   
83   for(s=1;s<255;s++){
84     /* P(s|1), prob of receiving s given 1 transmitted */
85     p1 = normal(((s-OFFSET+0.5)/amp - 1)/noise) -
86       normal(((s-OFFSET-0.5)/amp - 1)/noise);
87     
88     /* P(s|0), prob of receiving s given 0 transmitted */
89     p0 = normal(((s-OFFSET+0.5)/amp + 1)/noise) -
90       normal(((s-OFFSET-0.5)/amp + 1)/noise);
91     
92 #ifdef notdef
93     printf("P(%d|1) = %lg, P(%d|0) = %lg\n",s,p1,s,p0);
94 #endif
95     metrics[0][s] = gr_log2(2*p0/(p1+p0)) - bias;
96     metrics[1][s] = gr_log2(2*p1/(p1+p0)) - bias;
97   }
98   /* 255 is also a special value */
99   /* P(s|1) */
100   p1 = 1 - normal(((255-OFFSET-0.5)/amp - 1)/noise);
101   /* P(s|0) */
102   p0 = 1 - normal(((255-OFFSET-0.5)/amp + 1)/noise);
103   
104   metrics[0][255] = gr_log2(2*p0/(p1+p0)) - bias;
105   metrics[1][255] = gr_log2(2*p1/(p1+p0)) - bias;
106 #ifdef  notdef
107   /* The probability of a raw symbol error is the probability
108    * that a 1-bit would be received as a sample with value
109    * 0-128. This is the offset normal curve integrated from -Inf to 0.
110    */
111   printf("symbol Pe = %lg\n",normal(-1/noise));
112 #endif
113   for(bit=0;bit<2;bit++){
114     for(s=0;s<256;s++){
115       /* Scale and round to nearest integer */
116       mettab[bit][s] = floor(metrics[bit][s] * scale + 0.5);
117 #ifdef  notdef
118       printf("metrics[%d][%d] = %lg, mettab = %d\n",
119              bit,s,metrics[bit][s],mettab[bit][s]);
120 #endif
121     }
122   }
123 }