co_tpc.cpp

通信系统里面的信道编码中的乘积码合作编码visual c++程序

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


#define N 63
#define K 57
#define M (N-K+1)
#define Infolength K*K
#define pi 3.1416
double fm=0.02;
double fs=10;
double Mr=64;
int Nr=3696;
double dt=1/fs;

int H[N-K][N];
int G[K][N];
int a[N][N]={0};
int g[7]={1,1,0,0,1,1,1};
int r[N][N-K]={0};
int Alpha[64][6];
int k[Infolength];
int Info[K][K];
int Code[N][N];
int Num; 

double startsnr=6;
double endsnr=18;
double snrstep=1;

                                      //迭代次数
double Receive[N][N];
double alpha[8]={0.0,0.2,0.3,0.5,0.7,0.9,1.0,1.0};                           //重量因子
double beta[8]={0.2,0.4,0.6,0.8,1.0,1.0,1.0,1.0};                           //修正因子
int test[8][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,};

void input(int *message)
{
	int i;
	for(i=0;i<Infolength;i++)
		message[i]=rand()%2;
}

double UniformRand()
{
	double uRand=double(rand())/(double)RAND_MAX;
	return uRand;
}

double GaussNoise(double sigma)
{
	int num=12;
	double sum=0.0;
	double mean=0.0;
	double gaussRand;

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

double *Rayleigh(double fm,double Mr,double dt,int Nr)
{
	int i,j;
	double *tr;
	double *xr;
	double *yr;
	double *rr;
	tr = (double*)malloc(sizeof(double)*Nr);
	xr = (double*)malloc(sizeof(double)*Nr);
	yr = (double*)malloc(sizeof(double)*Nr);
	rr = (double*)malloc(sizeof(double)*Nr);

	for(i=0;i<Nr;i++)
	{
		tr[i]=dt*i;
	}
//	cout<<"a"<<endl;
	double c=sqrt(2/Mr);
//	cout<<c<<endl;
	double w=2*pi*fm;
	for(i=0;i<Nr;i++)
	{
		xr[i]=0;
		yr[i]=0;
		rr[i]=0;
	}
	for(i=0;i<Mr;i++)
	{
		double alpha=(2*pi*i-pi+(2*pi*UniformRand()-pi))/(4*Mr);
		double ph1=2*pi*UniformRand()-pi;
		double ph2=2*pi*UniformRand()-pi;
		for(j=0;j<Nr;j++)
		{
			xr[j]=xr[j]+c*cos(w*tr[j]*cos(alpha)+ph1);
			yr[j]=yr[j]+c*cos(w*tr[j]*sin(alpha)+ph2);
		}
	}
//	cout<<"b"<<endl;
	for(i=0;i<Nr;i++)
	{
		rr[i]=sqrt(xr[i]*xr[i]+yr[i]*yr[i])/sqrt(2);
	
	}
//	for(i=0;i<20;i++)
//	{
//		cout<<r[i*10]<<",  ";
//	}
//	cout<<endl;
//	cout<<"c"<<endl;
	
	
	free(tr);
	free(xr);
	free(yr);
	return rr;
	free(rr);
}

void Channel1(double sigma,int Code[N][N],double Receive[N][N])
{
	int i,j;
	double *z;
	z = (double*)malloc(sizeof(double)*Nr);
	z = Rayleigh(fm,Mr,dt,Nr);
	for(i=0;i<K;i++)
	{
		for(j=0;j<N;j++)
			Receive[i][j]=(double)(2*Code[i][j]-1)*z[i*j+j]/0.9+GaussNoise(sigma);
	}
//	cout<<"d"<<endl;
	free(z);
}

void Channel2(double sigma,int Code[N][N],double Receive[N][N])
{
	int i,j;
	double *z;
	z = (double*)malloc(sizeof(double)*Nr);
	z = Rayleigh(fm,Mr,dt,Nr);
	for(i=K;i<N;i++)
	{
		for(j=0;j<N;j++)
			Receive[i][j]=(double)(2*Code[i][j]-1)*z[i*j+j]/0.1+GaussNoise(sigma);
	}
//	cout<<"d"<<endl;
	free(z);
}
//计算平方欧氏距离
double DIS(double A3[][N],int A4[][N],int i,int j)
{
	int k;
	double dis=0;
	for(k=0;k<N;k++)
	{
		dis=(dis+((A3[i][k]-A4[j][k])*(A3[i][k]-A4[j][k])));
	}
	return dis;
}

double DIS_line(double A3[][N],int A4[][N],int i,int j)
{
	int k;
	double dis=0;
	for(k=0;k<N;k++)
	{
		dis=(dis+((A3[k][i]-A4[j][k])*(A3[k][i]-A4[j][k])));
	}
	return dis;
}

int Error_Statistic(int a[K][K],int b[K][K])
{
	int i,j,Enum;
	Enum=0;
	for(i=0;i<K;i++)
	{
		for(j=0;j<K;j++)
		{
			if(a[i][j]!=b[i][j])
				Enum++;
		}
	}
	return Enum;
}
//生成G,H矩阵
GH()
{
	for(int m=0;m<63;m++)
	{
		a[m][m]=1;
	}
	for(int i=0;i<63;i++)
	{
		for(int j=0;j<57;j++)
		{
			if(a[i][j]==0)
			{
			}
			else
			{
				for(int k=0;k<=6;k++)
				{
					a[i][j+k]=((a[i][j+k]+g[k])%2);
				}
			}
		}
		for(j=0;j<6;j++)
		{
			r[i][j]=a[i][j+57];
			Alpha[63-i][j]=a[i][j+57];
		}
	}
	for(i=K;i<N;i++)
	{
		H[i-K][i]=1;
	}
	for(i=0;i<K;i++)
	{
		G[i][i]=1;
	}
	for(i=0;i<57;i++)
	{
		for(int j=0;j<6;j++)
		{
			G[i][57+j]=r[i][j];
			H[j][i]=r[i][j];
		}
	}
	for(i=0;i<(M-1);i++)
	{
		for(int j=0;j<N;j++)
			cout<<H[i][j]<<",";
		cout<<endl;
	}
	cout<<endl; 
/*	
	FILE *fp;

	fp = fopen("H.txt", "a");
    if(fp == NULL)
    { 
		printf("\n open file error");
	    exit(0);
    } 
//输出H[N-K][N]	
	fprintf(fp, "\nH = [\n");
	for(i=0;i<(N-K);i++)
	{
		for(int j=0;j<N;j++)
			fprintf(fp,"%d,",H[i][j]);
		fprintf(fp,"\n");
	}
	fprintf(fp,"]\n");
//输出G[K][N]
	fprintf(fp, "\nG = [\n");
	for(i=0;i<K;i++)
	{
		for(int j=0;j<N;j++)
			fprintf(fp,"%d,",G[i][j]);
		fprintf(fp,"\n");
	}
	fprintf(fp,"]\n");
//输出Alpha[64][6]	
	fprintf(fp, "\nAlpha = [\n");
	for(i=0;i<64;i++)
	{
		for(int j=0;j<6;j++)
			fprintf(fp,"%d,",Alpha[i][j]);
		fprintf(fp,"\n");
	}
	fprintf(fp,"]\n");
	
	fclose(fp);
	*/
}

encode_bch_row(int i)
{
	int p,q;
	int code_bch[N]={0};
	for(p=0;p<N;p++)
	{
		for(q=0;q<K;q++)
		{
			code_bch[p]=(code_bch[p]+Info[i][q]*G[q][p])%2;
		}
		Code[i][p]=code_bch[p];
	}
	
}

encode_bch_line(int i)
{
	int p,q;
	int code_bch[N]={0};
	for(p=0;p<N;p++)
	{
		for(q=0;q<K;q++)
		{
			code_bch[p]=(code_bch[p]+Code[q][i]*G[q][p])%2;
		}
		Code[p][i]=code_bch[p];
	}
	
}

encode_TPC()
{
	int i,j;
	int t=0;
	input(k);
	for(i=0;i<K;i++)
	{
		for(j=0;j<K;j++)
		{
			Info[i][j]=k[t];
			t++;
		}
	}
//行编码
	for(i=0;i<K;i++)
	{
		encode_bch_row(i);
	}
//列编码
	for(i=0;i<N;i++)
	{
		encode_bch_line(i);
	}
}

decode_bch(int bch_code[8*N],int i)
{
	int p,q,s1,s2,t=0,error,sigma1;
	int S1[6]={0,0,0,0,0,0};
	int S2[6]={0,0,0,0,0,0};
	for(p=0;p<N;p++)
	{
		for(q=0;q<6;q++)
		{
			S1[q]=(S1[q]+Alpha[N-p][q]*bch_code[i*N+p])%2;
		}
	}
	for(p=0;p<6;p++)
	{
		t=(t+S1[p]);
	}
	if(t!=0)
	{
		t=0;
		for(p=0;p<N && t<6;p++)
		{
			t=0;
			for(q=0;q<6;q++)
			{
				if(S1[q]==Alpha[p+1][q])
				{
					t=t+1;
				}
			}
			s1=p;
		}
//		cout<<s1<<p<<endl;
		s2=(s1*2)%N;
//		cout<<s2<<endl;
		sigma1=-(s2-s1);
		error=(sigma1-1+N)%N;
//		cout<<error<<endl;
		bch_code[i*N+error]=(bch_code[i*N+error]+1)%2;
	}
	else
	{
//		cout<<i<<"行无错"<<endl;
	}
}

chase_decode_row()
{
	int i,j,k;
	for(i=0;i<K;i++)
	{
		int p1,p2,p3;
		int *test_code;
		test_code=(int*)malloc(sizeof(int)*(8*N));
//		int optional_code[8][N];
//		int *optional_code;
//		optional_code=(int*)malloc(sizeof(int)*(8*N));
		if(fabs(Receive[i][0])<fabs(Receive[i][1]))
		{
			if(fabs(Receive[i][1])<fabs(Receive[i][2]))
				p1=0,p2=1,p3=2;
			else
			{
				if(fabs(Receive[i][2])<fabs(Receive[i][0]))
					p1=2,p2=0,p3=1;
				else
					p1=0,p2=2,p3=1;
			}
		}
		else
		{
			if(fabs(Receive[i][0])<fabs(Receive[i][2]))
				p1=1,p2=0,p3=2;
			else
			{
				if(fabs(Receive[i][2])<fabs(Receive[i][1]))
					p1=2,p2=1,p3=0;
				else
					p1=1,p2=2,p3=0;
			}
		}	
	
		for(j=3;j<N;j++)
		{
			if(fabs(Receive[i][j])<fabs(Receive[i][p1]))
			{
				p3=p2;
				p2=p1;
				p1=j;
			}
			else
			{
				if(fabs(Receive[i][j])<fabs(Receive[i][p2]))
				{
					p3=p2;
					p2=j;
				}
				else
				{
					if(fabs(Receive[i][j])<fabs(Receive[i][p3]))
					{
						p3=j;
					}
				}
			}
		}
//		cout<<p1<<"---"<<p2<<"---"<<p3<<endl;
//		cout<<Code[i][p1]<<"---"<<Code[i][p2]<<"---"<<Code[i][p3]<<endl;
//		cout<<Receive[i][p1]<<"---"<<Receive[i][p2]<<"---"<<Receive[i][p3]<<endl;
		for(j=0;j<8;j++)
		{
			for(k=0;k<N;k++)
			{
				if(Receive[i][k]>=0)
				{
					test_code[j*N+k]=1;
				}
				else
				{
					test_code[j*N+k]=0;
				}
			}
		}
		for(j=0;j<8;j++)
		{
			test_code[j*N+p1]=test[j][0];
			test_code[j*N+p2]=test[j][1];
			test_code[j*N+p3]=test[j][2];
		}