📄 rm31_chase.c
字号:
// ------------------------------------------------------------------------
// File: rm31_chase.c
//
// Simulation of RM(3,1), equivalent to the (8,4,4) extended Hamming code.
// Soft-decision decoding performed by Chase type-II algorithm.
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only.
// If this program is used to perform simulations whose results are
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or
// liable for damage or loss that may be caused by the use of this program.
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#define MAX_RANDOM LONG_MAX // Maximum value of random()
#define RATE 0.5 // Coding rate = 4/8
#define INIT_SNR 3.0 // Initial value of Eb/N0
#define FINAL_SNR 8.0 // Final value of Eb/N0
#define SNR_INCREMENT 1.0 // Increment in Eb/N0
#define NUMSIM 1000000 // Number of simulations (one per 4 bits)
#define SWAP(a,b) itemp=(a); (a)=(b); (b)=itemp;
#define n 8
#define k 4
#define nk 4 // n-k
int L;
// Generator matric in systematic form
// Information positions: 1,2,3 and 4 (exchanged 4 and 5 in original)
//
int G[k][n] = { 1,0,0,0,1,1,1,0,
0,1,0,0,1,1,0,1,
0,0,1,0,1,0,1,1,
0,0,0,1,0,1,1,1 };
// Original:
// int G[k][n] = { 1,0,0,1,0,1,1,0,
// 0,1,0,1,0,1,0,1,
// 0,0,1,1,0,0,1,1,
// 0,0,0,0,1,1,1,1 };
int wh[16] = { 0, 1, 1, 2, 1, 2, 2, 3, /* Hamming weight function: */
1, 2, 2, 3, 2, 3, 3, 4 }; /* wh[i] = weight of i */
double sim, block_error;
double ber;
double amp;
double seed;
int error;
int data[k], codeword[n];
int data_int;
float snr;
float transmited[n];
float received[n];
int hard[n];
int decoded[n];
float reli[n];
float reli_ord[n+1]; // Indexed from 1 .. n
int perm[n+1]; // Indexed from 1 .. n
void initialize(void);
void awgn(void);
void encode(void);
void bpsk_dec(void);
void indexx(int ns, float arr[], int indx[]);
void dec2bin(int test, int err[]);
float metric(int vector[]);
void gen_err(int count[], int patt[]);
main()
{
int i,j;
unsigned long value, test;
int err[n], patt[n];
double max, temp;
int count[n];
int got_one;
printf("Enter the value of L: "); scanf("%d", &L);
// Compute 2^L - 1
value = 1;
for (i=0; i<L; i++) value = value << 1;
value -= 1;
printf("2^L-1 = %ld\n", value);
snr = INIT_SNR;
// *********** S I M U L A T I O N L O O P **************
while ( snr < (FINAL_SNR+0.001) )
{
initialize();
while (sim < NUMSIM)
{
// ----------- Random information
for (i=0; i<k; i++)
data[i] = (random()>>10) & 0x01;
/* convert data[] to integer for error computation purposes */
data_int = 0;
for (i=0; i<k; i++)
data_int = (data_int << 1) ^ data[i];
// ----------- RM encode
encode();
// ----------- BPSK mapping
for (i=0; i<n; i++)
if (codeword[i]) transmited[i] = -1.0;
else transmited[i] = 1.0;
// ----------- AWGN channel
awgn();
// ----------- Symbol-by-symbol decision --> h
bpsk_dec();
// ----------- Reliabilities of received symbols
for (i=0; i<n; i++)
reli[i] = fabs(received[i]);
for (i=0; i<n; i++)
reli_ord[i+1] = reli[i];
// Order reliabilities in ascending order
// ----------- Obtain permutation. NOTE: perm[1]..perm[n]
indexx(n, reli_ord, perm);
#ifdef DEBUG
printf("prm = ");
for (i=0; i<n; i++) printf("%4d ",perm[i+1]);
printf("\n");
#endif
test = 0;
for (i=0; i<n; i++) patt[i] = hard[i];
max = -DBL_MAX; // For metric computation
got_one = 0; // Codeword found indication
while (test <= value)
{
#ifdef DEBUG
printf("test = %d pattern = ",test);
for (i=0; i<n;i++) printf("%4d ",patt[i]);
printf(" check = %d\n",check(patt));
#endif
// ----------- Check that (h+e) is a codeword
if (check(patt))
{
// ----------- Compute likelihood, compare & choose best
got_one = 1;
#ifdef DEBUG
printf(" metric = %f\n", metric(patt));
#endif
if ((temp=metric(patt)) > max)
{
max = temp;
for (i=0; i<n; i++) decoded[i] = patt[i];
}
}
// ----------- Increase counter over least L reliable positions
test += 1; // (Could do GRAY COUNTING here)
dec2bin(test,count);
// ----------- Generate error pattern, e
gen_err(count,patt);
}
// ------------ Decode information
i = 0;
if (got_one) // If codeword was found
for (j=0; j<k; j++)
i = (i<<1) + decoded[j];
else // else output hard decision
for (j=0; j<k; j++)
i = (i<<1) + hard[j];
error = i ^ data_int;
/* if (error) block_error+=1.0; /* block error rate */
ber += (float) wh[error];
sim+=1.0;
}
printf("%f %8.0f %8.0f %13.8e\n", snr, ber, (k*sim), (ber/(sim*k)));
fflush(stdout);
snr += SNR_INCREMENT;
}
}
void encode()
{
int i,j;
for (j=0; j<n; j++)
{
codeword[j] = 0;
for (i=0; i<k; i++)
codeword[j] ^= ( data[i] * G[i][j] ) & 0x01;
}
#ifdef DEBUG
printf("\n-------\ncod = ");
for (j=0; j<n; j++)
printf("%4d ",codeword[j]);
printf("\n");
#endif
}
void bpsk_dec()
{
int j;
for (j=0; j<n; j++)
{
if (received[j] <= 0.0)
hard[j] = 1;
else
hard[j] = 0;
}
#ifdef DEBUG
printf("rec = ");
for (j=0; j<n; j++)
printf("%4.2lf ",received[j]);
printf("\n");
printf("hrd = ");
for (j=0; j<n; j++)
printf("%4d ",hard[j]);
printf("\n");
#endif
}
void awgn()
{
double u1,u2,s,noise,randmum;
int i;
for (i=0; i<n; 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] = transmited[i] + noise/amp;
#ifdef NO_NOISE
received[i] = transmited[i];
#endif
}
}
void indexx(int ns, float arr[], int indx[])
{
// Indexes an array arr[1..n]: Outputs indx[1..n] such that
// arr[indx[j]] is in ascending order. Neither n nor arr are changed.
//
{
int i, indxt, ir=ns, itemp, j,ks,l=1;
int jstack = 0, istack[50];
float a;
for (i=1; i <= ns; i++)
indx[i] = i;
for (;;)
{
if ( (ir-l) < 7) {
for (j=l+1; j <= ir; j++)
{
indxt = indx[j];
a = arr[indxt];
for (i=j-1; i >= 1; i--)
{
if (arr[indx[i]] <= a) break;
indx[i+1] = indx[i];
}
indx[i+1] = indxt;
}
if (jstack == 0) break;
ir = istack[jstack--];
l = istack[jstack--];
}
else {
ks = (l+ir) >> 1;
SWAP(indx[ks],indx[l+1]);
if (arr[indx[l]]>arr[indx[ir]]) {
SWAP(indx[l],indx[ir]);
}
if (arr[indx[l+1]]>arr[indx[ir]]) {
SWAP(indx[l+1],indx[ir]);
}
if (arr[indx[l]]>arr[indx[l+1]]) {
SWAP(indx[l],indx[l+1]);
}
i = l+1;
j = ir;
indxt = indx[l+1];
a = arr[indxt];
for (;;)
{
do i++; while (arr[indx[i]] < a);
do j--; while (arr[indx[j]] > a);
if (j < i) break;
SWAP(indx[i],indx[j]);
}
indx[l+1] = indx[j];
indx[j] = indxt;
jstack += 2;
if (jstack > 50) { printf("\nError in indxx\n"); exit(0); }
if ((ir-i+1) >= j-1) {
istack[jstack] = ir;
istack[jstack-1] = i;
ir = j-1;
}
else {
istack[jstack] = j-1;
istack[jstack-1] = l;
l = i;
}
}
}
}
}
void initialize()
{
time(&seed);
srandom(seed);
amp = sqrt(2.0*RATE*pow(10.0,(snr/10.0)));
block_error = 0.0;
ber = 0.0;
sim = 0.0;
}
int check(int patt[])
{
// Compute syndromes
int i,j;
int syn[nk];
for (i=0; i<nk; i++)
{
syn[i] = 0;
for (j=0; j<k; j++)
syn[i] ^= ( patt[j] * G[j][i+k] ) & 0x01;
syn[i] ^= patt[k+i];
if (syn[i])
return(0); // If syndrome nonzero test fails
}
return(1);
}
void dec2bin(int test, int count[])
{
int i;
for (i=0; i<L; i++)
{
count[i] = test%2;
test /=2;
}
}
float metric(int vector[])
{
int i;
float temp;
temp = 0.0;
for (i=0; i<n; i++)
temp += -(2.0*(float)vector[i]+1.0)*received[i];
return(temp);
}
void gen_err(int count[], int patt[])
{
// Generate error pattern from vector representing the combination of
// least reliable positions
int i;
for (i=0; i<n; i++) patt[i] = hard[i];
for (i=0; i<L; i++)
if (count[i])
patt[perm[i+1]-1] ^= 1;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -