📄 main_4_2_0.c
字号:
/* GLOBECOM PAPER SCHEME # 1
* File: main_4_2_0.c
* Date: August 17, 1996
* ----------------------------------------------------------------------
* Simulation of block coded 8PSK modulation using BCH codes of length 64
* To achieve graceful degradation over a satellite (AWGN) channel
* (for the 5th Mini-Conference in Communication Theory, GLOBECOM'96)
* ----------------------------------------------------------------------
*
* C_1 : (64,18,22) BCH code
* C_2 : (64,45, 8) BCH code
* C_3 : (64,63, 2) code
*
* Ordered statistics soft-decision decoding of the BCH codes, from Marc
* Fossorier routines.
*
* Order-4 soft-decision of C1 : BCH(64,18,22)
* Order-2 soft-decision of C2 : BCH(64,45, 8)
* Order-0 soft-decision of C3 : BCH(64,63, 2)
*/
// ------------------------------------------------------------------------
// 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 "def3.h"
/* #define NO_NOISE /* */
/* #define PRINT_DETAILS /* */
/* #define DEBUG /* */
#define WATCHIT /* */
#define n 64 /* Code length */
#define k1 18 /* Dimension first code */
#define k2 45 /* Dimension second code */
#define k3 63 /* Dimension third code */
#define RATE 1.96875 /* R = (18+51+63)/64 */
#define MAXRAND LONG_MAX /* for random number generation */
#define NUMSIM 1000 /* no. simulations per 100 codewords */
#define INIT 12.0
#define LIMIT 13.0
#define INC 0.5
/* Constant used in preprocessing for third stage */
#define SQR2 0.7071067812 /* sqrt(2)/2 */
main()
{
unsigned iloop, inner;
unsigned int codeword1[n], codeword2[n], codeword3[n];
unsigned int mess_one, mess_three;
unsigned int message1[k1], message2[k2], message3[k3];
unsigned int seed;
unsigned int word[n], estimate1[n], estimate2[n], estimate3[n];
unsigned int result, ytest;
int error1, error2, error3;
int G1[k1][n]; /* Generator matrix of BCH (64,18,22) */
int G2[k2][n]; /* Generator matrix of BCH (64,45, 8) */
int G3[k3][n]; /* Generator matrix of BCH (64,63, 2) */
double rx[n], ry[n], rx3[n], ry3[n];
double amp;
double pi = M_PI;
double error_one =0.0, error_two = 0.0, error_three = 0.0;
double q[2][8];
double rndm, u1, u2, s, x1, x2, aux, br;
double snr, pe, pe1, pe2, pe3;
register int i,j;
int fflush();
long int count_even = 0, count_odd = 0;
/* Variables for soft-decision based on ordered statistics. */
double out_D[N]; /* Decoded word: {-1,1} --> {0,1} */
double out_D2[N]; /* Decoded word: {-1,1} --> {0,1} */
double out_D3[N]; /* Decoded word: {-1,1} --> {0,1} */
void order4_dec(); /* Order-4 decoding of C1 : BCH(64,18,22) */
void order2_dec2(); /* Order-2 decoding of C2 : BCH(64,45, 8) */
void order0_dec3(); /* Order-0 decoding of C3 : BCH(64,63, 2) */
/* data file pointers */
FILE *fp1, *fp2, *fp3;
/* READ GENERATOR MATRIX OF BCH(64,18,22) CODE */
fp1 = fopen("generator_641822.data","r");
printf("\nGenerator matrix of BCH(64,18,22):\n");
for (i=0; i<k1; i++)
{
for (j=0; j<n; j++)
{
fscanf(fp1,"%1d", &G1[i][j]);
printf("%1d",G1[i][j]);
}
printf("\n");
}
fclose(fp1);
printf("\n");
/* READ GENERATOR MATRIX OF BCH(64,45, 8) CODE */
fp1 = fopen("generator_644508.data","r");
printf("\nGenerator matrix of BCH(64,45, 8):\n");
for (i=0; i<k2; i++)
{
for (j=0; j<n; j++)
{
fscanf(fp1,"%1d", &G2[i][j]);
printf("%1d",G2[i][j]);
}
printf("\n");
}
fclose(fp1);
printf("\n");
/* READ GENERATOR MATRIX OF BCH(64,63, 2) CODE */
fp1 = fopen("generator_646302.data","r");
printf("\nGenerator matrix of BCH(64,63, 2):\n");
for (i=0; i<k3; i++)
{
for (j=0; j<n; j++)
{
fscanf(fp1,"%1d", &G3[i][j]);
printf("%1d",G3[i][j]);
}
printf("\n");
}
fclose(fp1);
printf("\n");
fp1 = fopen("stage1-410-new.data","w");
fp2 = fopen("stage2-410-new.data","w");
fp3 = fopen("stage3-410-new.data","w");
snr = INIT;
/* Random seed from current time */
time(&seed);
srandom(seed);
do {
amp = sqrt( 2.0 * RATE * pow(10.0, (snr/10.0)) );
error_one = error_two = error_three = 0.0;
/* Compute q(0) and q(1), the in-phase and quadrature components
of the 8-PSK signal constellation --> BLOCK PARTITIONING (BP)
Label: (b1,b2,b3)
b1: Selects right/left half plane
b2: Selects upper/lower half plane
b3: Selects signal in a quadrant
|
100 o | o 000
|
101 o | o 001
|
-----------------------
|
111 o | o 011
|
110 o | o 010
|
*/
q[0][0]=amp*cos(3.0*pi/8.0);
q[1][0]=amp*sin(3.0*pi/8.0); /* 000 3pi/8 */
q[0][1]=amp*cos(pi/8.0);
q[1][1]=amp*sin(pi/8.0); /* 001 pi/8 */
q[0][2]=amp*cos(3.0*pi/8.0);
q[1][2]=-amp*sin(3.0*pi/8.0); /* 010 -3pi/8 */
q[0][3]=amp*cos(pi/8.0);
q[1][3]=-amp*sin(pi/8.0); /* 011 - pi/8 */
q[0][4]=-amp*cos(3.0*pi/8.0);
q[1][4]=amp*sin(3.0*pi/8.0); /* 100 5pi/8 */
q[0][5]=-amp*cos(pi/8.0);
q[1][5]=amp*sin(pi/8.0); /* 101 7pi/8 */
q[0][6]=-amp*cos(3.0*pi/8.0);
q[1][6]=-amp*sin(3.0*pi/8.0); /* 110 -5pi/8 */
q[0][7]=-amp*cos(pi/8.0);
q[1][7]=-amp*sin(pi/8.0); /* 111 -7pi/8 */
for (iloop=0; iloop<NUMSIM; iloop++)
{ /* begin big simulation loop */
error1 = error2 = error3 = 0;
for (inner=0; inner<100; inner++)
{ /* begin inner simulation loop */
/*
********************************************************************
ENCODING
********************************************************************
*/
/* --------- First-level code --------- */
for (j=0; j<k1; j++)
/* random k1-bit message */
codeword1[j] = message1[j] = (random()>>10)&01;
/* BCH (64,18,22) codeword */
for (j=k1; j<n; j++)
{
codeword1[j] = 0;
for (i=0; i<k1; i++)
codeword1[j] ^= (message1[i] & G1[i][j]);
}
/* --------- Second-level code -------- */
for (j=0; j<k2; j++)
/* random k2-bit message */
codeword2[j] = message2[j] = (random()>>10)&01;
/* BCH (64,45,8) codeword */
for (j=k2; j<n; j++)
{
codeword2[j] = 0;
for (i=0; i<k2; i++)
codeword2[j] ^= (message2[i] & G2[i][j]);
}
/* --------- Third-level code --------- */
for (j=0; j<k3; j++)
{
/* random k3-bit message */
codeword3[j] = message3[j] = (random()>>10)&01;
if (codeword3[j])
count_odd++;
else
count_even++;
}
/* BCH (64,63,2) codeword */
for (j=k3; j<n; j++)
{
codeword3[j] = 0;
for (i=0; i<k3; i++)
codeword3[j] ^= (message3[i] & G3[i][j]);
if (codeword3[j])
count_odd++;
else
count_even++;
}
/* --- 8-PSK Signal Mapping by Block Partitioning --- */
for (j=0; j<n; j++)
word[j] = 4.0*codeword1[j]+2.0*codeword2[j]+codeword3[j];
#ifdef PRINT_DETAILS
printf("\n");
for (j=0; j<n; j++)
printf("%d",word[j]);
printf("\n");
#endif
/*
********************************************************************
MODULATION + CHANNEL
********************************************************************
*/
for(j=0; j<n; j++)
{
#ifdef RAYLEIGH
do
{
rndm = (double)(random())/MAXRAND;
u1 = rndm * 2 - 1.0;
rndm = (double)(random())/MAXRAND;
u2 = rndm * 2 - 1.0;
s = u1 * u1 + u2 * u2;
} while( s >= 1 );
x1 = u1 * sqrt( (-2.0*log(s))/s );
x2 = u2 * sqrt( (-2.0*log(s))/s );
br = sqrt(x1*x1 + x2*x2); /* Rayleigh fading component */
#else
br = 1.0;
do
{
rndm = (double)(random())/MAXRAND;
u1 = rndm * 2 - 1.0;
rndm = (double)(random())/MAXRAND;
u2 = rndm * 2 - 1.0;
s = u1 * u1 + u2 * u2;
} while( s >= 1 );
x1 = u1 * sqrt( (-2.0*log(s))/s ); /* Gaussian component */
x2 = u2 * sqrt( (-2.0*log(s))/s ); /* Gaussian component */
/* Received (complex) signal */
aux = 2.0 * (double)word[j] * pi / 8.0 ;
rx[j] = br * q[0][word[j]] + x1 ; /* */
ry[j] = br * q[1][word[j]] + x2 ; /* */
#endif
#ifdef NO_NOISE
rx[j] = q[0][word[j]] ; /* */
ry[j] = q[1][word[j]] ; /* */
#endif
/* Save received symbols for third stage decoding */
rx3[j] = rx[j]; ry3[j] = ry[j];
#ifdef DEBUG
printf("%1d%1d%1d: rx[%d]=%lf, ry[%d]=%lf\n",
codeword1[j],codeword2[j],codeword3[j],
j,rx[j],j,ry[j]);
#endif
}
/*
********************************************************************
ORDERED STATISTICS DECODING
********************************************************************
*/
/* ------------- FIRST STAGE ------------ */
/* Order-4 soft-decision decoding of C1 : BCH(64,18,22) */
order4_dec(G1,rx,out_D);
/* Map {+1,-1} ---> {1,0} (2*bit - 1) */
for (j=0; j<n; j++)
if (out_D[j] == 1.0)
estimate1[j] = 1;
else
estimate1[j] = 0;
#ifdef PRINT_DETAILS
for (i=0;i<n;i++)
if (out_D[i]==1.0)
{ printf("1"); }
else
{ printf("0"); }
printf("\n");
#endif
/* ------------- SECOND STAGE ------------ */
/* Order-2 soft-decision decoding of C2 : BCH(64,45, 8) */
order2_dec2(G2,ry,out_D2);
/* {+1,-1} ---> {1,0} using (2*bit - 1) */
for (j=0; j<n; j++)
if (out_D2[j] == 1.0)
estimate2[j] = 1;
else
estimate2[j] = 0;
#ifdef PRINT_DETAILS
for (i=0;i<n;i++)
if (out_D2[i] == 1.0)
{ printf("1"); }
else
{ printf("0"); }
printf("\n");
#endif
/* ------------- THIRD STAGE ------------ */
/* Rotate and translate received signals according to the
decisions of the first two decoding stages. Projection
in the new x-axis */
for (i=0; i<n; i++)
if ((estimate1[i]==0) && (estimate2[i]==0))
/* 1st quadrant */
rx3[i] = - SQR2*(-ry3[i] + rx3[i]);
else
if ((estimate1[i]==0) && (estimate2[i]==1))
/* 4th quadrant */
rx3[i] = - SQR2*(ry3[i] + rx3[i]);
else
if ((estimate1[i]==1) && (estimate2[i]==1))
/* 3rd quadrant */
rx3[i] = SQR2*(-ry3[i] + rx3[i]);
else
if ((estimate1[i]==1) && (estimate2[i]==0))
/* 2nd quadrant */
rx3[i] = SQR2*(ry3[i] + rx3[i]);
/* Order-0 soft-decision decoding of C3 : BCH(64,63, 2) */
order0_dec3(G3,rx3,out_D3);
/* Map {+1,-1} ---> {1,0} (2*bit - 1) */
for (j=0; j<n; j++)
if (out_D3[j] == 1.0)
estimate3[j] = 1; /* */
else
estimate3[j] = 0; /* */
#ifdef PRINT_DETAILS
for (i=0;i<n;i++)
if (out_D3[i]==1.0)
{ printf("1"); }
else
{ printf("0"); }
printf("\n");
#endif
/* Count errors only in message positions */
for (j=0; j<k1; j++)
error1 += (codeword1[j]^estimate1[j]);
for (j=0; j<k2; j++)
error2 += (codeword2[j]^estimate2[j]);
for (j=0; j<k3; j++)
error3 += (codeword3[j]^estimate3[j]);
#ifdef DEB
for (j=0; j<n; j++) printf ("%d",codeword1[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",codeword2[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",codeword3[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",word[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",estimate1[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",estimate2[j]);
printf("\n");
for (j=0; j<n; j++) printf ("%d",estimate3[j]);
printf("\n---------------------------------------\n");
#endif
}
error_one += error1;
error_two += error2;
error_three += error3;
#ifdef DEBU
printf("%5d %10.0f %10.0f %10.0f\n", iloop, error_one,
error_two, error_three);
#endif
} /* end big simulation loop */
#ifdef WATCHIT
printf("Simulations (per signal sequence) = %g\n", (NUMSIM*100.0));
printf("Errors: %g %g %g %g\n", ( error_one + error_two + error_three ),
error_one, error_two, error_three);
#endif
/* Bit error rates */
pe1 = error_one/(NUMSIM*100.0*k1);
pe2 = error_two/(NUMSIM*100.0*k2);
pe3 = error_three/(NUMSIM*100.0*k3);
pe = ( error_one + error_two + error_three ) / (NUMSIM*100.0*(k1+k2+k3));
printf("%10.5f\t%10.7f\t%10.7f\t%10.7f\t%10.7f\n",
snr, pe, pe1, pe2, pe3);
fflush(stdout);
fprintf(fp1,"%10.5f\t%14.10f\n", snr, log10(pe1));
fprintf(fp2,"%10.5f\t%14.10f\n", snr, log10(pe2));
fprintf(fp3,"%10.5f\t%14.10f\n", snr, log10(pe3));
fflush(fp1);
fflush(fp2);
fflush(fp3);
snr += INC;
} while (snr <= LIMIT);
printf("counters, even and odd = %ld %ld\n",count_even, count_odd);
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -