spa.cpp

Implements layered decoding using the sum-product-algorithm. Codes are input as a .txt file in A-lis

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
//#include <time.h>
#include <ctime>

#define RATE 124.0/305.0 
#define PI M_PI
#define No 2
#define MAX_MN 25
#define MAX_MN_LENGTH 10000
#define K_SIZE 1260
#define MAX_BLOCK_ERRORS 100
#define MAX_BLOCKS 10000
#define HOURS 6
#define MAX_TIME HOURS*60.0*60.0
#define ITERATION_POINTS 10
#define SEED 546546513

double unif_rand(void);
double Vanessa(double,int);
double psi(double);
int sign(double);
double norm_rand(void);
double rayl1(void);
double phase(void);
double min(double,double);
double correction(double,double);
double Q(double,double);
double max(double, double);
int error_counter(int , int*,int*);

long int x = -1;


int main (int argc, char * const argv[]) 
{
	char out_file[25],in_file[25];
	FILE *file_pt;
	int v,m,n,max_n,max_m,data_length,temp,BSC,parity_ok,count,MIN,MULT_ITER,exclude,index,product;
	double eb_no,snr_number,sigma,epsilon,snr_coded,sum,sum_decision,previous_min;
	int pivot,iter,iterations,firstpass;
	unsigned long int bit_error_count,codeword_error_count;
	//double bit_error,codeword_error;
	double  start_snr,finish_snr,increment_snr,codeword_error_vector[ITERATION_POINTS],bit_error_vector[ITERATION_POINTS];
	unsigned long int codeword_error_count_vector[ITERATION_POINTS],bit_error_count_vector[ITERATION_POINTS],block_count,input_block;
	//int iter_value[ITERATION_POINTS];
	long double newly_received,previously_sent;
	int iter_value[ITERATION_POINTS] = {1,2,3,4,5,6,7,8,9,10};
	unsigned long int iter_average[ITERATION_POINTS];
	double average[ITERATION_POINTS];
	clock_t start,finish;
	double time;

	
    for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"seed"))
        {
            x = atoi(argv[v+1]);        //initializing seed;
            break;
        }
    }
	
	
    for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"snr"))
        {
            start_snr = atof(argv[v+1]);        
            finish_snr = atof(argv[v+2]);
            increment_snr = atof(argv[v+3]);
			break;
        }
    }
    
	
    for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"iter"))
        {
            iterations = strtoul(argv[v+1],NULL,10);
            break;
        }
    }        
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"in_file"))
        {
            
			strcpy(in_file,argv[v+1]);
			break;
        }
    }        
	
	
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"out_file"))
        {
            
			strcpy(out_file,argv[v+1]);
			break;
        }
    }        
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"BSC"))
        {
            BSC = atoi(argv[v+1]);
            break;
        }
	}        
	
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"block"))
        {
            input_block = strtoul(argv[v+1],NULL,10);
            break;
        }
    }        
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"MIN"))
        {
            MIN = atoi(argv[v+1]);
            break;
        }
	}        
	
	
	for(v=0;v<argc;v++)
    {
        if(!strcmp(argv[v],"MULT_ITER"))
        {
            MULT_ITER = atoi(argv[v+1]);
            break;
        }
	}        
	
	
    
	
    if(x == -1)                          //seed is not set
    {
        //printf("Enter seed number for random generator:  ");
        //scanf("%li",&x);
        x =  546465461; //time(0);
    }
	
	
    printf("Starting Seed: %li\n\n",x);
    /*Perform file I/O */
	
	
	file_pt = fopen(in_file,"r");
	
	if(file_pt == NULL)
	{
		printf("File not found \n");
		return 1;
	}
	
	fscanf(file_pt,"%d",&n);
      fscanf(file_pt,"%d",&m);
	data_length = -(m-n);
	
	//data_length = K_SIZE;
    fscanf(file_pt,"%d",&max_n);
    fscanf(file_pt,"%d",&max_m);
	
	//dynamically allocate memory
	
	int* n_length = new int [n];
	int* m_length = new int [m];
	int* parity = new int [m];
	int* actual_sent = new int [n];
	int* decoded = new int [n];
	long double* APP = new long double [n];
	long double* gamma = new long double [n];
	long double* received_signal = new long double [n];
	long double* rho = new long double [max_m];
	long double* output = new long double [max_m];
	long double* alpha = new long double [max_m];
	long double* beta = new long double [max_m];
	
	int** n_vector = new int* [n];
	int** m_vector = new int* [m];
	int** alpha_vi_fj = new int* [n];
	
	long double** Beta_vi_fj = new long double* [n];
	long double** vi_fj = new long double* [n];
	long double** fj_vi = new long double* [m];
	long double** lambda = new long double* [m];
	
	
	for(int k=0;k<n;k++)
	{
		n_vector[k] = new int [max_n];
		alpha_vi_fj[k] = new int [max_n];
		Beta_vi_fj[k] = new long double [max_n];
		vi_fj[k] = new long double [max_n];
	}
	
	for(int k=0;k<m;k++)
	{
		m_vector[k] = new int [max_m];
		fj_vi[k] = new long double [max_m];
		lambda[k] = new long double [max_m];
	}
	
	
	
	for(int i=0;i<n;i++)
	{
		fscanf(file_pt,"%d",&n_length[i]);
    }
	
	for(int i=0;i<m;i++)
	{
		fscanf(file_pt,"%d",&m_length[i]);
    }
	
	
	for(int i=0;i < n;i++)
	{
		for(int j = 0;j < n_length[i];j++)
		{
			fscanf(file_pt, "%d",&temp); 
			n_vector[i][j] = temp;	
		}
	}
    
	for(int i=0;i < m;i++)
	{
		for(int j = 0;j < m_length[i];j++)
		{
			
			
			fscanf(file_pt, "%d",&temp); 
			m_vector[i][j] = temp;
		}
	}
	
	
	
	printf("Eb_No_Start: \t %e \n", start_snr);
	printf("Eb_No_Finish: \t %e \n", finish_snr);
	printf("Eb_No_Increment: \t %e \n", increment_snr);
	
	if(MULT_ITER == 1)
	{
		printf("Running Multiple Iterations \n");
		iterations = iter_value[5];
	}
	else
	{
		printf("Iterations: \t %d \n", iterations);
	}
	
	printf("Blocks: \t %d \n", input_block);
	if(BSC == true)
	{
		printf("Performing Hard Decisions \n");
	}
	else
	{
		printf("Performing Soft Decisions \n");
	}
	
	switch(MIN)
	{
		case 0:
			printf("Not Performing Min-Sum Algorithm \n");
			break;
		case 1:
			printf("Performing Min-Sum Algorithm \n");
			break;
		case 2:
			printf("Performing Min-Sum Correction Algorithm \n");
			break;
		default:
			break;
	}
	
	
	
	printf("Input File: %s \n\n\n\n",in_file);
	
	start = clock();
	
	for(eb_no = start_snr; eb_no <= finish_snr; eb_no += increment_snr)
    {
		
		snr_number = pow(10.0,(double)eb_no/10.0);
		
		snr_coded = /*((double)data_length/(double)n)*/ RATE * snr_number;
		sigma = 1.0/sqrt(2.0*snr_coded);
		
		if (BSC == true)
		{
			epsilon = Vanessa(10.0 * log10(snr_coded),1000000);
		}
		
		bit_error_count = 0;
		codeword_error_count = 0;
		block_count = 0;
		for(int i=0;i<ITERATION_POINTS;i++)
		{
			codeword_error_count_vector[i] = 0;
			bit_error_count_vector[i] = 0;
			iter_average[i] = 0;
		}
		
		
		
		
		
while( (block_count < MAX_BLOCKS) || (eb_no >= 1.4 && codeword_error_count_vector[ITERATION_POINTS-1] <= 130))     
        {
			for(int k = 0; k < n; k++)
			{
				received_signal[k] =  1 + sigma * norm_rand();
				if(received_signal[k] > 0)
				{
					decoded[k] = 0;
				}
				else
				{
					decoded[k] = 1;
				}
				
				
				actual_sent[k] = 0;
				
				
				
				
				if(BSC == true)
				{   
					APP[k] =  pow(-1.0,actual_sent[k]) * log( (1.0-epsilon)/(epsilon) );
				}
				else
				{
					APP[k] = 2.0 * received_signal[k] / pow(sigma,2.0);
				}    
			}
			
				
				for(int i = 0;i < n ; i++)
				{                                                 //initialize vi -> fj = 0
						for(int k = 0; k < n_length[i]; k++)
						{
							vi_fj[i][k] = 0.0;
						}
						
				}

			
			
			for(iter = 0; iter <= iter_value[ITERATION_POINTS-1]; iter++)
			{
				
				for(int i=0;i < m;i++)
				{
					for(int k=0; k <m_length[i];k++)
					{
						sum = 0.0;
						
						for(int z=0;z<n_length[m_vector[i][k]-1];z++)
						{
							 if(n_vector[m_vector[i][k]-1][z]-1 != i)
							 {
								sum += P[m_vector[i][k]-1][z];										
							 }
						}
						
						Beta[k] = sum + APP[m_vector[i][k]-1];
						alpha[k] = sign(Beta[k]);
						Beta[k] = fabs(Beta[k]); 
					}
				
					for(int k=0;k<m_length[i];k++)
					{
						sum = 0.0;
						product = 1;
						
						int(int z=0;z<m_length[i];z++)
						{
							if(k != z)
							{
								sum += psi(Beta[z]);
								product = product * alpha[z];
							}
						}
						
						Q[m_vector[i][k]-1] = psi(sum) * product;
						
					}
				
					for(int k=0;k<m_length[i];k++)
					{
						for(int z =0; z< n_length[m_vector[i][k]-1];z++)
						{
							if(i == n_vector[m_vector[i][k]-1][z])
							{
								P[m_vector[i][k]-1][z] = Q[m_vector[i][k]-1];  
							}
						}
					}
				
				
				}
				
				for(int i=0;i<n;i++)
				{
					sum = 0;
					for(int k =0;k<n_length[i];k++)
					{
						sum += P[i][k];
					}
					
					sum_decision = sum + APP[i];
				}
				
								
							
									
					if( sum > 1e50 || sum < -1e50)
					{
						printf("Possible error \n");
					}
					
					
					if(sum_decision > 0)
					{
						decoded[i] = 0;
					}
					else
					{    
						decoded[i] = 1;
					}
					
				}
				
				//Perform Parity Check
				
				parity_ok = true;
				
				for(int i = 0;i<m;i++)
				{
					parity[i] = 0;    
					for(int k=0;k<m_length[i];k++)
					{
						pivot = m_vector[i][k] - 1; 
						parity[i] = (parity[i] + decoded[pivot]) % 2;    
					}