spa.cpp

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

					if(parity[i] == true)
					{ 
						parity_ok = false;
						break;
					}
				}
	
				for(int i = 0;i < ITERATION_POINTS;i++)
				{
					if(iter  == iter_value[i])
					{
						count = error_counter(n,decoded,actual_sent);
						if(count != 0)
						{
							codeword_error_count_vector[i]++;
						}
						bit_error_count_vector[i] += count;
						iter_average[i] += iter;
					}
				}
					
																		
																											
																																			
				if(parity_ok == true && iter != 0)           //don't continue iterations because parity checks
				{
					iter++;                 //need to update iteration count since breaking
					break;
				}
				
				
				
			}	
			
			
			
			count = error_counter(n,decoded,actual_sent);
			
			iter--;      //decrementing for the increment that the for loop causes
			
			for(int i =0;i<ITERATION_POINTS;i++)                     //account for missed iterations
			{
				if(i == 0)
				{
					if(iter<iter_value[i])
					{
						for(int c = i; c<ITERATION_POINTS;c++)
						{
							if(count != 0)
							{
								codeword_error_count_vector[c]++;
							}
							bit_error_count_vector[c] += count;
							iter_average[c] += iter; 
						}
					}
				}	
				else if(iter >= iter_value[i-1] && iter < iter_value[i])
				{
					for(int c=i; c<ITERATION_POINTS;c++)
					{
						if(count != 0)
						{
							codeword_error_count_vector[c]++;
						}
						bit_error_count_vector[c] += count;
						iter_average[c] += iter;
					}
				}
					
			}
			


			finish = clock();
			time = (double(finish)-double(start))/CLOCKS_PER_SEC;
			
			if(time > MAX_TIME)
			{	
				
				printf("Time: %f \n\n\n\n",time);
				break;
			}




			block_count++;	
		}
		
		
		for(int c=0;c<ITERATION_POINTS;c++)
		{
			
			bit_error_vector[c] = (double)bit_error_count_vector[c]/(double)(n/*data_length*/*block_count);
			codeword_error_vector[c] = (double)codeword_error_count_vector[c]/(double)block_count;
			average[c] = (double)iter_average[c]/(double)block_count;
		}
		
		if(BSC == true)
		{
			
			for(int c=0;c<ITERATION_POINTS;c++)
			{
				printf("%f \t %e \t %e \t %e \n",eb_no,epsilon,bit_error_vector[c],codeword_error_vector[c]);
			}
		}
		else
		{
			
			printf("iter      ");
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %d               ",iter_value[c]);
			}
			printf("\n \n");

			printf("iter_aver");

                        for(int c =0;c<ITERATION_POINTS;c++)
                        {
                                printf("\t %e ",average[c]);
                        }

			
			printf("\n");
			
			printf("bit_e     ");
			
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %e ",bit_error_vector[c]);
			}
			printf("\n");
			
			printf("codeword_e");
			
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %e ",codeword_error_vector[c]);
			}
			printf("\n \n \n");
			
			
			printf("iter      ");
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %d   ",iter_value[c]);
			}
			printf("\n \n");
			
			printf("bit_c     ");
			
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %d ",bit_error_count_vector[c]);
			}
			printf("\n");
			
			printf("codeword_c");
			
			for(int c =0;c<ITERATION_POINTS;c++)
			{
				printf("\t %d ",codeword_error_count_vector[c]);
			}
			printf("\n \n \n");
			
			
			printf("Block_Count: %d \n",block_count);
			
			
			
			
			
			//printf("iter \t bit_e \t \t \t codeword_e \n");
			
			/*		for(int c=0;c<ITERATION_POINTS;c++)
			 {
			 printf("%d \t \t  %e \t %e \n",iter_value[c],bit_error_vector[c],codeword_error_vector[c]);
			 }
			 
			 printf("\n \n");
			 printf("iter \t bit_c \t codeword_c\n");
			 
			 for(int c=0;c<ITERATION_POINTS;c++)
			 {
			 printf("%d \t \t %u \t %u \n",iter_value[c], bit_error_count_vector[c],codeword_error_count_vector[c]);
			 }
			 */		
		}
		
	}
	
	
	
	
	
	
	fclose(file_pt);
	
	delete [] n_length;
	delete [] m_length;
	delete [] parity;
	delete [] actual_sent;
	delete [] decoded;
	delete [] APP;
	delete [] gamma;
	delete [] received_signal;
	delete [] rho;
	delete [] alpha;
	delete [] beta;
	delete [] output;
	
	
	
	for(int k=0;k<n;k++)
	{
		delete n_vector[k];
		delete alpha_vi_fj[k];
		delete Beta_vi_fj[k];
		delete vi_fj[k];
	}
	
	for(int k=0;k<m;k++)
	{
		delete fj_vi[k];
		delete lambda[k];
		delete m_vector[k];
	}
	
	
	
	delete [] n_vector;
	delete [] m_vector;
	delete [] alpha_vi_fj;
	delete [] Beta_vi_fj;
	delete [] vi_fj;
	delete [] fj_vi;
	delete [] lambda;
	
	return 0;
}


//------------------------------------------------------------------------------------------------------------------------------------------
//Returns a uniform random variable from 0...1
//------------------------------------------------------------------------------------------------------------------------------------------

double unif_rand(void)
{
	
    long int IM;
	
    IM = x/127773;
    x = 16807 * (x - IM * 127773)-IM *2836;
	
    if ( x < 0)
    {
        x = x + 2147483647;
    }
	
    return (x/(double)2147483647);
	//return (double)rand()/(double)RAND_MAX;
}

//---------------------------------------------------------------------------------------------------------------------------------
//Generates a Uniform random variable between (0,2*PI)
//---------------------------------------------------------------------------------------------------------------------------------
double phase(void)
{
    return(2.0*PI*unif_rand());
}
//----------------------------------------------------------------------------------------------------------------------------------
//Generates a Rayleigh random variable
//------------------------------------------------------------------------------------------------------------------------------------------
double rayl1(void)
{
    return(sqrt(-2.0*log(1.0-unif_rand())));
}
//----------------------------------------------------------------------------------------------------------------------------------
//Generates a Guassian R.V.  N(0,1)
//---------------------------------------------------------------------------------------------------------------------------------
double norm_rand(void)
{
    double y,Nrayl,Nphase;
    Nrayl=rayl1();
    Nphase=phase();
    y=Nrayl*cos(Nphase);
    return(y);
}


double Vanessa(double snr,int iterations)
{
	
	double snr_number,sigma,error,number,signal;
	int received, data;
	
	snr_number = pow(10,snr/10.0);  
    error = 0;
	
	
    sigma = 1/sqrt(2*snr_number);                                                     //standard deviation
	
	
    for(int j=0;j < iterations;j++)
	{
        number = unif_rand();                                                        //Gemerate Data
        
		if (number > .5)
		{
			data = 0;
        }
		else
		{
			data = 1;
        }
		
        
        if (data == 0)
		{
            signal = 1 + sigma * norm_rand();
        }
		else
        {    
			signal = -1 + sigma * norm_rand();
        }
        
        
        if (signal > 0)
		{
            received = 0;
        }
		else 
		{ 
			received = 1;
			
		}
		
        
        
        if (received != data)
		{
            error++;
        }
        
    }    
    
    
    return (double)error/(double)iterations;
	
}	

double psi(double x1)
{
	
	long double psi;
	long double c1;
	
	c1 = exp(x1);
	
	if(c1 == HUGE_VAL)
	{
		psi = 0;
	}
	else if(c1 == 1)
	{
		psi = 100000000000.0;
	}
	else
	{
		
		psi = log( (c1 +1.0)/(c1-1.0));
	}
	
	
	return psi;
}

int sign(double num)
{
	int test;
	
	if(num > 0)
	{
		test = 1;
	}
	else
	{
		test = -1;
	}
	
	return test;
}



double min(double x1,double x2)
{
	if(x1 < x2)
	{
		return x1;
	}
	else
	{
		return x2;
	}
	
	
}


double correction(double x1, double x2)
{
	if(fabs(x1+x2) < 2 && fabs(x1-x2) > 2 *fabs(x1+x2))
	{
		return .5;
	}
	else if(fabs(x1-x2) < 2 && fabs(x1+x2) > 2 * fabs(x1-x2))
	{
		return -.5;
	}
	else
	{
		return 0;
	}
}

int error_counter(int length , int* vector,int* actual)
{
	int r = 0;
	
	for(int q=0;q<length;q++)
	{      
		if( vector[q] != actual[q] )
		{
			r++;
		}
	}
	return r;
}

double max(double x1,double x2)
{
	if(x1 < x2)
	{
		return x2;
	}
	else
	{
		return x1;
	}
}


double Q(double x2,double y2)
{
	double value;
	
	value = max(x2,y2) + max(5.0/8.0 - fabs(x2 - y2)/4.0,0) - max(x2 + y2,0) - max(5.0/8.0 - fabs(x2 + y2)/4.0,0); 
	
	return value;
}