vitwgn_wimax2.cpp

针对WImax的实体层作设计,里面包括卷积编码,以及相对应的Viterbi解码器

#include<iostream.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define size 2048
#define pi   3.14159265
#define Nstate 4   //states
#define Nlink 2    //routines
#define Maxtry 1   //the number of try
#define Nfft 1024  //FFT points
#define Ncpc 4  // cyclic prefix for 16QAM
#define M 16    //M-ary alphabet
#define T2 Nfft*Ncpc //Ncpcs
#define T T2/2
#define M2 M/2

int sitlv=Ncpc/2;
int u[8]={1,3,-1,-3,5,-5,7,-7};
int v[8]={1,3,-1,-3,5,-5,7,-7};
int Aset[8][2][4]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
				   0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
				   0,4,1,5,2,6,3,7,8,12,9,13,10,14,11,15,
				   0,8,1,9,2,10,3,11,4,12,5,13,6,14,7,15};
double wr[Nfft],wi[Nfft]; //noise
int symr[M],symi[M];  //modulated symbols
int nn[T][2],kk[T][2];
double dt[T][2][2];
double Anoise=0;  //noise amplitude

int sb[Nstate][2]; //previous state binary value
int cnt[Nstate]={-1,-1,-1,-1}; //count
int ss[Nstate][2];
int ns[Nstate][2]; //next state
int ps[Nstate][2]; //previus state
int as[Nstate][2]; //data in
int cc[Nstate][2][2],a[T],c[T][2];
int ss1[Nstate];
int s,s0;

double xr[Nfft],xi[Nfft];//transmitted symbols
double Eb=0;
int cbit=0;
int ek[16]={0};

void fsm(int bit);
void dec2bin();
void interleaver();
void SDtable();
void modulation();
void VitWGN();
void WGN();
int Vit(int a);
int minimum(int x, int y);

void fft(int sign);

void main()
{
    srand(time(NULL));

	fsm(0);
	fsm(1);
	interleaver();
	SDtable();
	modulation();
	VitWGN();

	getchar();
}

//fsm
void fsm(int bit)
{
	int i,j;
	for (i = 0; i <Nstate; i++) {
		cnt[i]=-1;
	}
	dec2bin(); //10->2

	switch(bit){
	case 0:
	   for ( i = 0; i < Nstate; i++) {
		  ss[i][5]=bit;
		  ss[i][4]=sb[i][5];
		  ss[i][3]=sb[i][4];
		  ss[i][2]=sb[i][3];
		  ss[i][1]=sb[i][2];
		  ss[i][0]=sb[i][1];
		  //bin2dec
		  ss1[i]=ss[i][5];
		  for (j = 4; j >=0; j--) {
		  ss1[i]=2*ss1[i]+ss[i][j];
		  }
		  ns[i][bit]=ss1[i];
		  cnt[ss1[i]]+=1;
		  ps[ss1[i]][cnt[ss1[i]]]=i;
		  as[ss1[i]][cnt[ss1[i]]]=bit;
		  cc[i][0][bit]=fmod(bit+sb[i][0]+sb[i][1]+sb[i][2]+sb[i][5],2);
		  cc[i][1][bit]=fmod(bit+sb[i][0]+sb[i][1]+sb[i][2]+sb[i][4]+sb[i][5],2);
	   }
	   break;

	case 1:
	   for (i = 0; i < Nstate; i++) {
		  ss[i][5]=bit;
		  ss[i][4]=sb[i][5];
		  ss[i][3]=sb[i][4];
		  ss[i][2]=sb[i][3];
		  ss[i][1]=sb[i][2];
		  ss[i][0]=sb[i][1];
		  ss1[i]=ss[i][5];
		  //bin2dec
		  for (j = 4; j >=0; j--) {
		  ss1[i]=2*ss1[i]+ss[i][j];
		  }
		  ns[i][bit]=ss1[i];
		  cnt[ss1[i]]+=1;
		  ps[ss1[i]][cnt[ss1[i]]]=i;
		  as[ss1[i]][cnt[ss1[i]]]=bit;
		  cc[i][0][bit]=fmod(bit+sb[i][0]+sb[i][1]+sb[i][2]+sb[i][5],2);
		  cc[i][1][bit]=fmod(bit+sb[i][0]+sb[i][1]+sb[i][2]+sb[i][4]+sb[i][5],2);
	   }
	   break;
	}
}

void dec2bin()
{
	int i,n,d,k;

	for (k = 0; k < Nstate; k++) {
		  if (k==0) {
			 for (n = 0; n <= 5; n++) {
				sb[k][n]=0;
			 }
		  }
		  else {
			 d=k;
			 for ( i=0; d!=0; i++ ) {
				 sb[k][i] = d % 2;
				 d = d / 2 ;
			 }

			 for (n = i+1; n <= 5; n++) {
				 sb[k][n]=0;
			 }
		   }
	   }
}

void interleaver()
{
	int d=16;  //interleaver parameter
	int em[T2]={0};
	int dm[T2]={0};
	int dk[T2]={0};

	for (int k = 0; k < T2; k++) {
		em[k]=(T2/d)*fmod(k,d)+floor(k/d);
		ek[k]=sitlv*floor(em[k]/sitlv)+fmod( em[k]+T2-floor(d*em[k]/T2),sitlv)+1;
		dm[k]=sitlv*floor(k/sitlv)+fmod( k+floor(d*k/T2),sitlv );
		dk[k]=d*dm[k]-(T2-1)*floor(d*dm[k]/T2)+1;
	}
}

void SDtable()
{
	int mv[M2][3]={0,0,0,0,0,1,0,1,0,0,1,1,1,0,0,1,0,1,1,1,0,1,1,1};
	int t,m,k,mk,pos;

	//Tables k and n
	for (t = 0; t < T; t++) {
	   for (m = 0; m < 2; m++) {
		   pos=ek[cbit]-1;
		   nn[t][m]=(int)(pos/Ncpc)+1;
		   kk[t][m]=fmod(pos,Ncpc);
		   cbit=cbit+1;
	   }
	}

	
}

void modulation()
{

	int j=0;
	double sum=0;

	switch (M) {

	   case 2:   //Bpsk
		 for (int m = -1; m <= 1; m+=2) {
			 symr[j]=m;
			 symi[j]=0;
			 j++;
		 }
		 break;

	   case 4:  //Qpsk,4QAM
		  for ( int m = 0; m <= 1; m++ ) {
			 for ( int n = 0; n <= 1; n++) {
				sum+=pow(u[m],2)+pow(v[n],2);
				symr[j]=u[m];
				symi[j]=v[n];
				j++;
			 }
		  }
		  Eb=sum/2*4;
		  break;

	   case 16:  //16QAM
		  for ( int m = 0; m <= 3; m++ ) {
			 for ( int n = 0; n <= 3; n++) {
				  sum+=pow(u[m],2)+pow(v[n],2);
				  symr[j]=u[m];
				  symi[j]=v[n];
				  j++;
			 }
		  }
		  Eb=sum/(4*16);
		  break;

	   case 64:  //64QAM
		  for ( int m = 0; m <= 7; m++ ) {
			 for ( int n = 0; n <= 7; n++) {
				  sum+=pow(u[m],2)+pow(v[n],2);
				  symr[j]=u[m];
				  symi[j]=v[n];
				  j++;
			 }
		  }
		  Eb=sum/(6*64);
		  break;
	}
}

void VitWGN()
{
	int z[T];
	int SNRperbit[10]={0};
	int np=0;
	double err;
	double Ex;
	double Ew;
	int ci[T][2];
	double rr[Nfft],ri[Nfft];
	double Pb[T];
	int t,k,m,n,c1,c2;

	for (double p = 5; p<= 10; p+=5) {
        cout<<p<<"\n";
		err=0;
		Ex=0;
		Ew=0;
		Anoise=sqrt(0.5*Eb/pow(10,p/10));

		for (int iter = 0; iter < Maxtry; iter++) {

		   //1/2 convolutional encoder
		   for (t = 0; t < T; t++) {
			   a[t]=0+(rand()%2);
		   }
		   for (k = 0; k < Nstate/2; k++) {
			   a[k]=a[k+T-Nstate/2]; //Randomizer
		   }
		   s0=( a[T-1])*2+a[T-2];
		   s=s0;

		   //encoder
		   for (t = 0; t < T; t++) {
			   c[t][0]=cc[s][0][a[t]];
			   c[t][1]=cc[s][1][a[t]];
			   s=ns[s][a[t]];
		   }

		   //interleaver
		   for (t = 0; t < T; t++) {
			   ci[t][0]=c[ek[t]][0];
			   ci[t][1]=c[ek[t]][1];
		   }

		   //modulation
		   int value[Nfft];
		   for (int i = 0; i < Nfft; i++) {
			   value[i]=0+(rand()%M);
		   }
		   for (n = 0; n < Nfft; n++) {
			   xr[n]=symr[value[n]];
			   xi[n]=symi[value[n]];
		   }
		   
		   //IFFT
		   fft(1);

		   //AWGN
		   WGN();
		   for (int n = 0; n < Nfft; n++) {
			   rr[n]=xr[n]+wr[n];
			   ri[n]=xi[n]+wi[n];
		   }

		   for (int n = 0; n < Nfft; n++) {
			   xr[n]=rr[n];
			   xi[n]=ri[n];
		   }
		   //FFT
		   fft(-1);

		   //demodulation
		   double D[M][Nfft];
		   for (n = 0; n < Nfft; n++) {
			  for (m = 0; m < M; m++) {
				  D[m][n]=sqrt(pow( (xr[n]-u[m]),2 )+pow( (xi[n]-v[m]),2 ));
			  }
		   }

		   double tt;
		   for (t = 0; t < T; t++) {
			  for (c1 = 0; c1 < 2; c1++) {
				 for (c2 = 0; c2 < 2; c2++) {
					 tt=0;
					 for (m = 0; m < M2; m++) {
						tt+=D[ Aset[m][0][kk[t][0]] ][ nn[t][0] ]+D[ Aset[m][1][kk[t][1]] ][ nn[t][1] ];
					 }
					 dt[t][c1][c2]=tt;
				 }
			  }
		   }

           //Viterbi Decoder
		   for (t = 0; t < T; t++) {
			   z[t]=Vit(t);
		   }

		   //error correction result
		  for (t = 0; t < T; t++) {
			 if (a[t]!=z[t]) {
				err+=(a[t]+z[t]);
			 }
		  }
		  //cout<<err<<"\n";
		  
		}
		Pb[np]=err/(T*Maxtry);
		cout<<Pb[np]<<"\n";
		np+=1;
		cout<<"\n";
	}

}

void WGN() //noise generator
{
		   const static int q = 15;
		   const static float c1 = (1 << q) - 1;
		   const static float c2 = ((int)(c1/3 )) + 1;
		   const static float c3 = 1.f / c1;
		   float random = 0.f;
		   float noise = 0.f;

		   for (int n = 0; n < Nfft; n++) {
			   random = ((float)rand() / (float)(RAND_MAX + 1));
			   noise = (2.f * ((random * c2) + (random * c2) + (random * c2)) - 3.f * (c2 - 1.f)) * c3;
			   wr[n]=noise;
			   wr[n]*=Anoise;
		   }
		   for (int n = 0; n < Nfft; n++) {
			   random = ((float)rand() / (float)(RAND_MAX + 1));
			   noise = (2.f * ((random * c2) + (random * c2) + (random * c2)) - 3.f * (c2 - 1.f)) * c3;
			   wi[n]=noise;
			   wi[n]*=Anoise;
		   }
}

int Vit(int a)
{
	int ac[Nstate], dist[Nlink];
	int acc[Nstate]={1000,1000,1000,1000};
	int pk[T][4]={0};
	int ct[4];
	int b[16];
	int sum=0;
	int t,i,j,k;

	acc[s0]=0; //ideal tailbiting
	for (t=0; t<T; t++) {
		for (i=0; i<4; i++) {
			for (k=0; k<Nlink; k++) {
				ct[0]=cc[0][ps[i][k]][as[i][k]];
				ct[1]=cc[1][ps[i][k]][as[i][k]];
				dist[k]=acc[ ps[i][k] ];
			}
			ac[i]=minimum(dist[0],dist[1]);
			if (ac[i]==dist[0]) {  //distinguish the previous state
				pk[t][i]=0;
			}
			else {
				pk[t][i]=1;
			}
		}
		for (j = 0; j <4; j++) {
             acc[j]=ac[j];
		}
	}

	s=s0; //ideal tailbiting
	for (t = T-1; t>=0; t--) {
		k=pk[t][s];
		b[t]=as[s][k];
		s=ps[s][k];
	}
	return b[a];
}

int minimum(int x, int y)
{
	int min;

	if (x<min) {
	   min=x;
	}
	if (y<min) {
	   min=y;
	}
	
	return min;
}

void fft(int sign)  //Computing FFT
{
	 long m, n, irem, l, le, le1, k, ip, i, j;
	 double ur, ui, wr, wi, tr, ti, temp;

	 j=0;  //Bit Reversal

	 for(i=0; i<Nfft-1; i++){
		if(i<j){
		   tr=xr[i];
		   ti=xi[i];
		   xr[i]=xr[j];
		   xi[i]=xi[j];
		   xr[j]=tr;
		   xi[j]=ti;
		   k=Nfft/2;
		   while(k<=j){
			  j=j-k;
			  k=k/2;
		  }
		}
		else{
		   k=Nfft/2;
		   while(k<=j){
			  j=j-k;
			  k=k/2;
		   }
		}
		j=j+k;
	 }

	 m=0;irem=Nfft; //Butterfly