turbo_arq.cpp

该程序是用Ｃ语言编写的

//this is the simulation program for RCPC hybrid ARQ
//384 bits of a frame and 1/4 parent code rate
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define PI 3.1415926
#define RAND_MAX 32767
void crcencode(void);        //crc encoder
void ptencoding(void);       //Turbo encoder
void punct(void);            //puncture processing
void depunct(void);          //depuncturing processing
void uncoded(void);          //uncoded
void decoding(void);         //tubor encoder
void noise(int);             //channel
void cdecode(void);          //crc decoder
void selective(void);        //arq process
void sranint(void);        //s-random interleaving
//double sranint(void);        //s-random interleaving
double rnd(void);            //random data
double norndl(void);         //gaussian random data

double sigma,sum=0,temp;
int w=1,l,n=368,m=384,nm=1356,row,repeat,rpt,mm=8,pun,pun1;
const int bps=34000,tail=384;
int cnum,t,nn=8,NN=34;    //parameters for rayleigh fading

//variables for tubor encoding
int tt1,tt2,tt3,tt4,vk,yy1[2500],yy2[2500],yy3[2500],yy4[2500],yy[6500],yin[6500];
int yz[2500],sg[6500],punc[6500],y[2500];

//variables for tubor decoding
int decision[2500];
float dpunc1[2500],dpunc2[2500],dpunc3[2500],dpunc4[2500],dpunc5[2500];
float intv1[2500],intv2[2500],Le21[2500],deleav[2500],dt[6500],dt1[6500],Lre21[2500],Lre12[2500];
float sigma1[2500][20],sigma2[2500][20],Lre1[2500],final[2500];
float dsigma1[2500][20],dsigma2[2500][20],dalpha[2500][20],dbeta[2500][20];
float pred[200][2500],alpha[2500][20],beta[2500][20],Lre2[2500],Lre[2500];

int zz[100][2000],si,ei,dd,vc,done,aa;   //crc and arq parameters
int p[]={1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1};  //crc polynomial

//variables & parameters for s-random interleaving
int sy[2500],sz[2500],sz1[2500],syy[2500],msb[2500],lsb[2500],oput[2500];
int dc,ii,tt,idex,idex1,index[2500],count,reverse[2500],kk=17;
float indat[2500],sintv[2500],slre[2500],smlre[2500],sran[2500],slre21[2500];

//table index for the s-random interleaving when the frame size is blow 512
int invalue[]={5,15,5,15,1,9,9,15,13,15,3,15,5,13,15,9,3,1,3,15,1,13,1,9,15,11,3,15,5};
int er,error,er1,error1,celer,celer1;   //variables for error counting



void main()
{
	int i=0,j,ri,retn;
	float tput[900],derror[900],sn,loop,retnm[900],celerr[900],effny[900];
    FILE *out1;
	out1=fopen("gtnfat.txt","w");
	loop=100;  //the number of frame simulated
//	tail=384;     //Frame size including tail bits
	l=10;         //first dB
    aa=1;         //counting block transmitting
	srand(time(NULL));

	//SNR range
	while(l<160)
	{
         //Initialization
		ri=0;         //first frame
		er=0;        //error rate without error correction
		celer1=0;    //frame error
		er1=0;       //bit error
		cnum=0;      //symbol duration
		t=0;         //counting frame duration
		retn=0;      //retransmission number

        //Frames at each SNR point

		while(ri<loop)
		{
			si=0;      //start index of window size
			ei=0;      //end index of windoe size
			done=0;
			vc=aa;     //initial value for counting the number fo frame
			dd=2;      //first puncturing pattern

			crcencode();  //crc encoding
			ptencoding();   //turbo encoding

			//ARQ process
		//	while(!done)

			if(!done)
			{
				punct();   //puncturing
				noise(l);  //channel
				depunct(); //depuncturing

				switch(dd)
				{
					//select depunctured parity bits
				case 1:{
					uncoded();
					break;
					   }
				case 2:{
					decoding();  //1/2 coding
					break;
					   }
				case 3:{
					decoding();   //1/3 coding
					break;
					   }
				case 4:{
					decoding();   //1/4 coding
					break;
					   }

				}
//////////////////
	        //	er=error+er;

             //   printf("er=:%6d\n",er);

		        cdecode();       //crc decoding	
		        selective();      //arq
			}

			er1=er1+error1;   //bit error rate counting
			celer1=celer1+celer;   //frame error rate counting
			retn=retn+vc;    //throughput counting
			ri++;
           
		}


sn=10*log10(1/((0.002+0.001*l)*(0.002+0.001*l)));   //SNR
//  Residual error rates

        derror[i]=er1/(loop*m);   //bit error rates
		celerr[i]=celer1/loop;    //frame error rates
//throughput
		retnm[i]=retn/(loop);     //average retransmission number
		tput[i]=0.958/retnm[i];   //throughput considered a overhead by crc and tail bits
//(384/16000)=0.024,line efficiency:distance 1500km-->0.005
		//distance 3000km-->0.01
		effny[i]=0.024/((0.024+0.01)*retnm[i]);
		printf("%f %6.5f %8.6f %6.5f %5.4f %5.4f\n",sn,tput[i],derror[i],celerr[i],retnm[i],effny[i]);
		fprintf(out1,"%f %6.5f %8.6f %6.5f %5.4f %5.4f\n",sn,tput[i],derror[i],celerr[i],retnm[i],effny[i]);
		if(tput[i]<0.1)
			break;
//SNR range inteval adjusted
		if(sn>25)
			l=l+6;
		if(sn<25&&sn>20)
			l=l+12;
		if(sn<20&&sn>15)
			l=l+24;
		if(sn<15&&sn>10)
			l=l+36;
		if(sn<10&&sn>7)
			l=l+48;
		if(sn<7)
			l=l+84;
/////////////
		i++;
//////////
	}
	fclose(out1);
	
}

void crcencode()
{
	int i,j,c;
	float x[384];
	int s[17];
	for(j=0;j<w;j++)//w=1
	{
		//random input data generation
		for(i=0;i<n;i++)
		{
            x[i]=rnd();
			if(x[i]>0.5)
				x[i]=1;
			if(x[i]<0.5)
				x[i]=0;

		}
		for(i=n;i<384;i++)
            x[i]=0;
		for(i=0;i<384;i++)
			y[i]=x[i];
//crc operation
		for(c=1;c<n+1;c++)
		{
			if(y[0]==1)
			{
				for(i=0;i<17;i++)
					s[i]=y[i]^p[i];
				for(i=0;i<17;i++)
					y[i]=s[i+1];
				for(i=17;i<384;i++)
					y[i-1]=y[i];
			}
			else if(y[0]==0)
			{
				for(i=0;i<17;i++)
					y[i]=y[i+1];
				for(i=17;i<384;i++)
                    y[i-1]=y[i];
			}
		}
//syndrome added
		for(i=0;i<16;i++)
			y[i+n]=y[i];
		for(i=0;i<n;i++)
			y[i]=x[i];
	}
}


//Turbo encoder
void ptencoding()
{
	int i,j,k;
	int x1[384];
	//encoder without interleaving
	//output encoded from the first recursive symtematic encoder
	for(pun=0;pun<m;pun+=mm)
	{
		tt1=0;
		tt2=0;
		tt3=0;
		if(y[0+pun]==0)
			yy1[0+pun]=0;
		if(y[0+pun]==1)
			yy1[0+pun]=1;
		tt1=yy1[0+pun];
		for(i=1;i<mm;i++)
		{
			vk=(tt3^tt2^y[i+pun]);
			yy1[i+pun]=vk^tt1^tt3;
			tt3=tt2;
			tt2=tt1;
			tt1=vk;
		}
	}
	//output encoded from the second sysmatic encoder
	for(pun=0;pun<m;pun+=mm)
	{
		tt1=0;
		tt2=0;
		tt3=0;
        if(y[0+pun]==0)
			yy2[0+pun]=0;
		if(y[0+pun]==1)
			yy2[0+pun]=1;
		tt1=yy2[0+pun];
		for(i=1;i<mm;i++)
		{
				vk=(tt3^tt2^y[i+pun]);
			yy2[i+pun]=vk^tt1^tt2^tt3;
			tt3=tt2;
			tt2=tt1;
			tt1=vk;
		}
	}
	//s-random interleaving
	sranint();
	for(i=0;i<m;i++)
		x1[i]=y[index[i]];
	//encoder with interleaving
	//output from the first RSC encoder
	for(pun=0;pun<m;pun+=mm)
	{
		tt1=0;
		tt2=0;
		tt3=0;
		if(x1[0+pun]==0)
			yy3[0+pun]=0;
		if(x1[0+pun]==1)
			yy3[0+pun]=1;
		tt1=yy3[0+pun];
		for(i=1;i<mm;i++)
		{
			vk=(tt3^tt2^x1[i+pun]);
			yy3[i+pun]=vk^tt1^tt3;
			tt3=tt2;
			tt2=tt1;
			tt1=vk;
		}
	}
	//output encoded from the second RSC encoder
    for(pun=0;pun<m;pun+=mm)
	{
		tt1=0;
		tt2=0;
		tt3=0;
		if(x1[0+pun]==0)
			yy4[0+pun]=0;
		if(x1[0+pun]==1)
			yy4[0+pun]=1;
		tt1=yy4[0+pun];
		for(i=1;i<mm;i++)
		{
			vk=(tt3^tt2^x1[i+pun]);
			yy4[i+pun]=vk^tt1^tt2^tt3;
			tt3=tt2;
			tt2=tt1;
			tt1=vk;
		}
	}
	//encoding data multiplexed
	for(i=0;i<m;i++)
	{
		yy[i]=y[i];
		yy[i+m]=yy1[i];
        yy[i+2*m]=yy3[i];
		yy[i+3*m]=yy4[i];
	}
}


//puncturing processer
void punct()
{
	int i,j;
	if(dd==1)
	{
		//code rate 1
		for(i=0;i<m;i++)
		{
			yz[i]=yy[i];
		}
	}
	if(dd==2)
	{
		//code rate 1/2
		for(i=0;i<192;i++)
		{
			yz[i*2]=yy[i*2+m];
			yz[i*2+1]=yy[i*2+1+2*m];
		}
	}
	if(dd==3)
	{
		//code rate 1/3
		for(i=0;i<192;i++)
		{
			yz[i*2]=yy[i*2+2*m];
			yz[i*2+1]=yy[i*2+1+m];

		}
	}
	if(dd==4)
	{
		//code rate 1/4
		for(i=0;i<m;i++)
		{
			yz[i]=yy[i+m*3];

		}
	}
	for(i=0;i<tail;i++)
	{
		punc[i]=yz[i];
		//Transmitting BPSK signals to channel
		if(punc[i]==1)
			punc[i]=1;
		if(punc[i]==0)
			punc[i]==-1;
	}
}


//rayleigh fading channel
void noise(int l)
{
	int i,j,dc,jj;
	float noise[tail],z[tail],n1[bps],n2[bps],temp,r1[bps];
	float ry1[bps],ry2[bps],ry3[bps],ry[tail],r2[bps],dop,fm,k1,k2;

	fm=185.0/16000.0 ;  //Doppler frequency/data transmission rate
		for(i=0;i<tail;i++)
		{
			//gaussian noise generation
            temp=norndl();
            noise[i]=(0.02*temp)+(0.001*l);
			cnum=cnum+1;
		}

//choose memoryless or correlated rayleigh fading
//memoryless fading
		for(i=0;i<tail;i++)
		{
			n1[i]=rnd();
			n2[i]=rnd();
		}
		for(i=0;i<tail;i++)
		{
			ry1[i]=sqrt(-2*log(n1[i]))*cos(2*PI*n2[i]);
			ry2[i]=sqrt(-2*log(n1[i]))*sin(2*PI*n2[i]);
		}
        for(i=0;i<tail;i++)
			ry[i]=sqrt(ry1[i]*ry1[i]+ry2[i]*ry2[i])/1.249891;
//correclated rayleigh fading
//Jake simulator
		for(i=0;i<tail;i++)
		{
			r1[i+t]=0;
			r2[i+t]=0;

			for(jj=1;jj<nn;jj++)
			{
				dop=cos(2*PI*fm*(cos(2*PI*jj/NN))*(i+t));
				r1[i+t]=r1[i+t]+cos(PI*jj/nn)*dop;
                r2[i+t]=r2[i+t]+sin(PI*jj/nn)*dop;
			}
			ry1[i+t]=(2*r1[i+t]+(sqrt(2)*cos(2*PI*fm*(i+t))))/(sqrt(nn));
			ry2[i+t]=(2*r2[i+t])/(sqrt(nn+1));
//185HZ:0.843766,5.5HZ:0.885909,for normalization
			ry3[i+t]=(sqrt(ry1[i+t]*ry1[i+t]+ry2[i+t]*ry2[i+t]))/sqrt((2)*0.843766);


		}
//checking one period
		for(i=0;i<tail;i++)
		
			ry[i]=ry3[i+t];
			if(cnum<16000)
				t=t+tail;
			else if(cnum>16000)
			{
				cnum=0;
				t=0;
			}
			for(i=0;i<tail;i++)
			{
				z[i]=(ry[i]*punc[i])+noise[i];
			}
/*			for(i=0;i<tail;i++)
			{
				if(z[i]>0.0)
					dt1[i]=1;
				if(z[i]<0.0)
					dt1[i]=-1;
			}
//data transmitted through channel
			for(i=0;i<tail;i++)
				dt[i]=z[i];
//checking the bit error rate without error correction
			error=0;
			for(i=0;i<tail;i++)
			{
				if(dt1[i]!=punc[i])
					error=error+1;
			}
*/		
}

//depuncturing and uncoding processors
void depunct()
{
	int i,j;
	if(dd==1)
	{
		//code rate 1
	for(i=0;i<m;i++)
		dpunc1[i+1]=dt[i];
    }
	if(dd==2)
	{
		//code rate 1/2 
		for(i=0;i<192;i++)
		{
			dpunc2[i*2+1]=dt[2*i];
			dpunc2[2*i+2]=0;
			dpunc3[2*i+2]=dt[2*i+1];
			dpunc3[2*i+1]=0;
			dpunc4[i+1]=0;
			dpunc5[i+1]=0;
		}
	}
	if(dd==3)
	{
		//code rate 1/3
		for(i=0;i<192;i++)
		{
			dpunc3[2*i+1]=dt[2*i];
			dpunc2[2*i+2]=dt[2*i+1];
		}
	}
	if(dd==4)
	{
		//code rate 1/4
		for(i+0;i<m;i++)
            dpunc5[i+1]=dt[i];
	}
}

void uncoded()
{
	int i,j;
	for(j=si;j<w-ei;j++)
	{
		for(i=0;i<m;i++)
		{
			if(dpunc1[i+1]>0)
				zz[j][i]=1;
			if(dpunc1[i+1]<0)
				zz[j][i]=0;
		}
    }
    for(i=0;i<tail;i++)
	{
		if(dpunc1[i]>0.0)
			dt1[i]=1;
		if(dpunc1[i]<0.0)
			dt1[i]=-1;
	}
    error1=0;
    celer=0;
	for(i=0;i<tail;i++)
	{
		if(dt1[i]!=punc[i])
		  error1=error1+1;
	}
	if(error1!=0)
		celer++;		
}