sccc.cpp

这是一个学习SCCC的很好的c 的程序

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

#define Infolength      1024
#define CodeRate        0.5         
#define numState_con        4
#define numState_dpsk      2
#define numBranch       2
#define ERRORNUM     100
#define InterLength    2048


int CodeWords0[Infolength * 2];
int CodeWords[Infolength * 2];
int InfoBit[Infolength];
double Receive[Infolength * 2];
int DecInfo[Infolength];
double sigma;
int InterPattern[InterLength];

double  startsnr =1;
double	endsnr = 9;
double	snrstep=0.5;


struct TransitionBranch
{
	int checkBit;
	int next_state;
}trellis_con[numState_con][numBranch]={0,0,1,2,0,2,1,0,1,3,0,1,1,1,0,3};

struct Transition
{
	int next_state;
}trellis_dpsk[numState_dpsk][numBranch]={1,0,0,1};

void sccc_dpsk_decoder_bcjr(double *receive, double *priori, double *extrinsic);
void sccc_con_decoder_bcjr(double *systematic, double *parity_check, double *priori, double *extrinsic_sys, double *extrinsic_par);
void simulation();
/////////////////////////////////////////////////////////////////////////////////////////

void Input(int *message)
{
	int i;

	for (i = 0; i < Infolength; i++)
		message[i] = rand() % 2;

}

void con_encoder()
{
	int i; 
	int CheckBit[Infolength];
	
	int ss = trellis_con[0][InfoBit[0]].next_state;

    CheckBit[0] = trellis_con[0][InfoBit[0]].checkBit;
    for(i = 1; i < Infolength; i++)
	{
		CheckBit[i] = trellis_con[ss][InfoBit[i]].checkBit;
		ss = trellis_con[ss][InfoBit[i]].next_state;
	}
	for(i = 0; i < Infolength; i++)
	{
		CodeWords0[2 * i] = InfoBit[i];
		CodeWords0[2 * i + 1] = CheckBit[i];
	}
}

void dpsk_encoder()
{
	int i;

	CodeWords[0] = trellis_dpsk[1][CodeWords0[0]].next_state;
	for (i = 1; i < Infolength * 2; i++)
	{
		CodeWords[i] = trellis_dpsk[CodeWords[i - 1]][CodeWords0[i]].next_state;	
	}

}
///////////////////////////////////////////////////////////////////////////////////////
void GenInterPattern()
{
	int i, j, test[InterLength];

	for(i = 0; i< InterLength; i++)
		test[i] = 0;
	for(i = 0; i < InterLength; i++)
	{
		j = rand() % InterLength;
		while (test[j] == 1)	
			j = rand() % InterLength;
		test[j] = 1;
		InterPattern[i] = j;
	}
}


void drandinter(double *interinput)
{
	int i;
	double test[InterLength];

	for(i = 0;i < InterLength; i++)
		test[InterPattern[i]] = interinput[i];
	for(i = 0;i < InterLength; i++)
		interinput[i] = test[i];

}

void randinter(int *interinput)
{
	int i;
	int test[InterLength];

	for(i = 0;i < InterLength; i++)
		test[InterPattern[i]] = interinput[i];
	for(i = 0;i < InterLength; i++)
		interinput[i] = test[i];

}

void dranddeinter(double *interinput)
{
	int i;
	double test[InterLength];

	for(i = 0;i< InterLength; i++)
		test[i] = interinput[InterPattern[i]];
	for(i = 0;i < InterLength; i++)
		interinput[i] = test[i];
}
////////////////////////////////////////////////////////////////////////////////////
double UniformRand()
{
    double uRand = double(rand())/(double)RAND_MAX;	 
    return uRand;
}

double GaussNoise(double sigma)
{						// sigma -- standard deviation of Gaussian noise
    int num = 12 ;            // total number of random to generate gaussRand
    double sum = 0.0;         // sum of uniformly distributed random number
    double mean = 0.0;  // mean of Gaussian noise
	double gaussRand;

    // compute sum of 12 uniform random number
    for (int i = 0; i < num; i++)
        sum += UniformRand();
    gaussRand = sigma * (sum - 6.0) * sqrt((double)num / 12.0) + mean;
    return gaussRand;
}

void Channel(double sigma, int *CodeWords, double *Receive)   //BPSK
{
	int i;
	for( i = 0 ; i < Infolength * 2; i++)
	{		
		Receive[i] = (double)(2 * CodeWords[i] - 1) + GaussNoise(sigma);
	}	
}
//////////////////////////////////////////////////////////////////////////////////////
void sccc_decoder()
{
	int i, num;
	double *Ys, *Yp;         
	double *La, *Le, *Les, *Lep;
	double Lc;            
	int Maxiter = 2;


    La = (double*)malloc(sizeof(double) *(Infolength * 2));   //动态分配存储块
	Le = (double*)malloc(sizeof(double) *(Infolength * 2));
	Les = (double*)malloc(sizeof(double) *Infolength); 
	Lep = (double*)malloc(sizeof(double) *Infolength);
	Ys = (double*)malloc(sizeof(double) *(Infolength * 2));
	Yp = (double*)malloc(sizeof(double) *Infolength );
	
    Lc = 2 / sigma / sigma;  
	for (i = 0; i < Infolength * 2 ; i++)
	{
		Le[i] = 0;
	}
    
	for (num = 0; num < Maxiter; num++)
	{
		for (i = 0; i < Infolength * 2 ; i++)
		{
			Ys[i] = Receive[i];
		    La[i] = Le[i];
		}

	    sccc_dpsk_decoder_bcjr(Ys, La, Le);
        dranddeinter(Le);
		for(i = 0; i < Infolength; i++)
		{
		    Ys[i] = Le[2 * i];
		    Yp[i] = Le[2 * i + 1];
		    La[i] = 0;
		}
	    sccc_con_decoder_bcjr(Ys,Yp,La,Les,Lep);
		if (num != Maxiter - 1)
		{
			for (i = 0; i < Infolength; i++)
			{
				Le[2 * i] = Les[i];
			    Le[2 * i + 1] = Lep[i];
			}
			drandinter(Le);
		}
	}

    for (i = 0; i < Infolength; i++)
	{
		if (Les[i] + Ys[i] > 0)
			DecInfo[i] = 1;
		else
			DecInfo[i] = 0;
	}
   
	free(La);
	free(Le);
	free(Les);
	free(Lep);
	free(Ys);
    free(Yp);
	
}    
////////////////////////////////////////////////////////////////////////////////////////
void sccc_dpsk_decoder_bcjr(double *receive, double *priori, double *extrinsic)
{
	int i, j, k;
	int sf, sb, br;
	double temp, Lc;            //  Lc channal reliability value
	double *reg;          // 
	double **afa, **beta, **gama;  

    Lc = 2 / sigma / sigma;                         //   Lc = 4 / N0;
	reg = (double*)malloc(sizeof(double)*numState_dpsk); 
	afa = (double**)malloc(sizeof(double)*(Infolength * 2 + 1));              //Length of afa is  Infolength + 1
	for (i = 0; i < Infolength * 2 + 1; i++)
		afa[i] = (double*)malloc(sizeof(double)*numState_dpsk);            // each has numState states
	beta = (double**)malloc(sizeof(double)*(Infolength * 2 +1));
    for (i = 0; i < Infolength * 2 + 1; i++)
		beta[i] = (double*)malloc(sizeof(double)*numState_dpsk);
	gama = (double**)malloc(sizeof(double)*(Infolength * 2));                       //Length of gama is  Infolength 
	for (i = 0; i < Infolength * 2; i++)
		gama[i] = (double*)malloc(sizeof(double)*(numState_dpsk * numBranch));     //there are only numState * numBranch different states in the trillis graph each time 

	// initialization afa, beta
	afa[0][1] = 1;
	afa[0][0] = 0;
	beta[Infolength][1] = 0.5;
	beta[Infolength][0] = 0.5;
	
	// calculate afa
	
	for(k = 0; k < Infolength * 2 ; k++)
	{
		for (i = 0; i < numState_dpsk; i++)
			reg[i] = 0;
		temp = 0.0;
		for(sf = 0; sf < numState_dpsk; sf++)                
		{
			for(br = 0; br < numBranch; br++)          
			{
				sb = trellis_dpsk[sf][br].next_state;       
				j = sf * numBranch + br;                
				gama[k][j] = exp(0.5 * (  (2 * trellis_dpsk[sf][br].next_state - 1) * Lc * receive[k] 
					                    + (2 * br - 1) * priori[k]));
		        reg[sb] += afa[k][sf] * gama[k][j]; 
				temp += afa[k][sf] * gama[k][j];
			}
		}                                                                           
	     for(i = 0; i < numState_dpsk; i++)
		    afa[k + 1][i] = reg[i] / temp;
	}

    // calculate beta
	for (k = Infolength * 2 - 1; k > 0; k--)
	{
		for (i = 0; i < numState_dpsk; i++)
			reg[i] = 0;
		temp = 0.0;
		for (sf = 0; sf < numState_dpsk; sf++)
		{
			for (br = 0; br < numBranch; br++)
			{
				sb = trellis_dpsk[sf][br].next_state;
				j = sf * numBranch + br;
				reg[sf] += beta[k+1][sb] * gama[k][j];
				temp += afa[k][sf] * gama[k][j];
			}
		}
	    for(i = 0; i < numState_dpsk; i++)
		    beta[k][i] = reg[i] / temp;        
	}



    for (k = 0; k < Infolength * 2; k++)
	{
	  	for (i = 0; i < numBranch; i++)
		     reg[i] = 0;
     	for (sf = 0; sf < numState_dpsk; sf++)
		{
	    	for (br = 0; br < numBranch; br++)
			{
			    sb = trellis_dpsk[sf][br].next_state;
			    j = sf * numBranch + br;
		    	reg[br] += afa[k][sf] * gama[k][j] * beta[k + 1][sb];
			}
		}
		extrinsic[k] = log(reg[1] / reg[0]) - priori[k];
	}
		

	free(reg);
	for(i = 0; i < Infolength * 2; i++)
		free(gama[i]);
	free(gama);
	for(i =0 ; i < (Infolength * 2 + 1); i++)
		free(afa[i]);
	free(afa);

	for(i = 0; i < (Infolength * 2 + 1); i++)
		free(beta[i]);
	free(beta);
}

///////////////////////////////////////////////////////////////////////

void sccc_con_decoder_bcjr(double *systematic, double *parity_check, double *priori, double *extrinsic_sys, double *extrinsic_par)
{
	int i, j, k;
	int sf, sb, br, c;
	double temp, Lc, *reg, *reg1;                                                       //  Lc channal reliability value         
	double **afa, **beta, **gama; 
	
    Lc = 2 / sigma / sigma;                                           //   Lc = 4 / N0;
	reg = (double*)malloc(sizeof(double)*numState_con); 
	reg1 = (double*)malloc(sizeof(double)*numBranch);
	afa = (double**)malloc(sizeof(double)*(Infolength + 1));          //Length of afa is  Infolength + 1
	for (i = 0; i < Infolength + 1; i++)
		afa[i] = (double*)malloc(sizeof(double)*numState_con);            // each has numState states
	beta = (double**)malloc(sizeof(double)*(Infolength +1));
    for (i = 0; i < Infolength + 1; i++)
		beta[i] = (double*)malloc(sizeof(double)*numState_con);
	gama = (double**)malloc(sizeof(double)*Infolength);                       //Length of gama is  Infolength 
	for (i = 0; i < Infolength; i++)
		gama[i] = (double*)malloc(sizeof(double)*(numState_con * numBranch));     //there are only numState * numBranch different states in the trillis graph each time 


	// initialization afa, beta
	afa[0][0] = 1;
	beta[Infolength][0] = 1.0 / numState_con;
	for (i = 1; i < numState_con; i++)
	{
        afa[0][i] = 0;
		beta[Infolength][i] = 1.0 / numState_con;
	}


	// calculate afa
	for(k = 0; k < Infolength ; k++)
	{
		for (i = 0; i < numState_con; i++)
			reg[i] = 0;
		temp = 0.0;
		for(sf = 0; sf < numState_con; sf++)                
		{
			for(br = 0; br < numBranch; br++)           
			{
				sb = trellis_con[sf][br].next_state;       
				j = sf * numBranch + br;                
				c = 2 * trellis_con[sf][br].checkBit - 1; 
				gama[k][j] = exp(0.5 * ((2 * br - 1) * (priori[k] +  systematic[k])+
					                      parity_check[k] * c));                     
		        reg[sb] += afa[k][sf] * gama[k][j];  
				temp += afa[k][sf] * gama[k][j];
			}
		}                                                                          
	   for(i = 0; i < numState_con; i++)
		   afa[k + 1][i] = reg[i] / temp;
	}


    // calculate beta
	for (k = Infolength - 1; k > 0; k--)
	{
		for (i = 0; i < numState_con; i++)
			reg[i] = 0;
		temp = 0.0;
		for (sf = 0; sf < numState_con; sf++)
		{
			for (br = 0; br < numBranch; br++)
			{
				sb = trellis_con[sf][br].next_state;
				j = sf * numBranch + br;
				reg[sf] += beta[k+1][sb] * gama[k][j];
				temp += afa[k][sf] * gama[k][j];
			}
		}
	    for(i = 0; i < numState_con; i++)
		    beta[k][i] = reg[i] / temp;        
	}



    for (k = 0; k < Infolength; k++)
	{
		for (i = 0; i < numBranch; i++){
		     reg[i] = 0;
			 reg1[i] = 0;
		}

     	for (sf = 0; sf < numState_con; sf++)
		{
	    	for (br = 0; br < numBranch; br++)
			{
			    sb = trellis_con[sf][br].next_state;
			    j = sf * numBranch + br;
		    	reg[br] += afa[k][sf] * gama[k][j] * beta[k + 1][sb];
				reg1[trellis_con[sf][br].checkBit] += afa[k][sf] * gama[k][j] * beta[k + 1][sb];
			}
		}
	    extrinsic_sys[k] = log(reg[1] / reg[0])  -  systematic[k];
		extrinsic_par[k] = log(reg1[1] / reg1[0])  -  parity_check[k];
	}
		

	free(reg);
	free(reg1);
	for(i = 0; i < Infolength; i++)
		free(gama[i]);
	free(gama);
	for(i =0 ; i < (Infolength + 1); i++)
		free(afa[i]);
	free(afa);
	for(i = 0; i < (Infolength + 1); i++)
		free(beta[i]);
	free(beta);
}
///////////////////////////////////////////////////////////////////////////////////////
int Error_Statistic()
{
	int i, Enum;
   
	Enum = 0;
	for (i = 0; i < Infolength; i++)
		if (DecInfo[i] != InfoBit[i])
			Enum++;
    return Enum;
}
//////////////////////////////////////////////////////////////////////////////////////////
void main()
{
	simulation();
}

void simulation()
{
	double snr, ber = 0.0, EbN0_dB;
	int errorNum ;
	int frameNum ;

	FILE *fp;

	fp = fopen("result.txt", "a");
    if(fp == NULL)
    { 
		printf("\n open file error");
	    exit(0);
    } 	
    
    fprintf(fp, "\n snr_db = [");
    for(EbN0_dB = startsnr;EbN0_dB <= endsnr;EbN0_dB += snrstep)
        fprintf(fp,"%f   ",EbN0_dB);
    fprintf(fp,"]\nber = [");
    fclose(fp);

    GenInterPattern();
	for(snr = startsnr;snr <= endsnr;snr += snrstep)
	{
		srand((unsigned) time(NULL));
		sigma = sqrt(pow(10,-0.1*snr)); // rate = 0.5
        
		errorNum = 0;
		frameNum = 0;
		while( errorNum < ERRORNUM)
		{
	        Input(InfoBit);
			con_encoder();
			randinter(CodeWords0);
	        dpsk_encoder();
            Channel(sigma, CodeWords, Receive);
			sccc_decoder();
	       	errorNum += Error_Statistic();
	      	frameNum++;
		}
		printf("%d", frameNum);
		printf("\n");
		ber = (double)errorNum / ( frameNum * Infolength);
	    printf("snr:%3.2lf, ber:%e \n", snr, ber);
    	printf("\n");
    	printf("\n");


		fp = fopen("result.txt","a");
        if(fp == NULL)
		{
			printf("\n open file error");
	        exit(0);
		}	
        fprintf(fp, "%e ", ber);
	    fclose(fp);

	}
	fp = fopen("result.txt","a");
    if(fp == NULL)
	{ 
		 printf("\n open file error");
	     exit(0);
	}
    fprintf(fp,"]\n");
    fclose(fp);


}
////////////////////////////////////////////////////////////////////////////////////////////////