📄 st.cpp
字号:
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include "Global.h"
#include "ST.h"
extern IO IO;
extern Random Rnd;
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
void CNM::Init( int length1, int length2 )
{
Length1 = length1;
Length2 = length2;
Rate = (float) length2 / (float) length1;
}
void CNM::SetChipInfo( float set, float *input )
{
int i;
for( i=Length1; i<Length2; i++ )
input[i] = set;
}
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
void CI::Init( int packSize )
{
PackSize = packSize;
MemoryAllocation( );
ChipInterleaverSet( );
}
void CI::MemoryAllocation( )
{
ShaffleRule = new int [PackSize];
Buffer = new NumericType [PackSize];
}
void CI::End( )
{
delete [] ShaffleRule;
delete [] Buffer;
}
void CI::ChipInterleaverSet( )
{
int i,j;
for ( j=0; j<NT; j++ )
{
for( i=0; i<PackSize; i++ )
ShaffleRule[i] = i;
Scramble( &ShaffleRule[0], PackSize, _ScrambleNum+NT*NT ); //??????????
for( i=0; i<PackSize; i++ )
Interleaver_ST[j][i]=ShaffleRule[i];
}
}
void CI::Interleaver( int *input )
{
int i,j;
for( i=0; i<PackSize; i++ ) Buffer[i].i = input[i];
// for( i=0; i<PackSize; i++ ) input[ShaffleRule[i]] = Buffer[i].i; //???????????
for ( j=0; j<NT; j++ )
for( i=0; i<PackSize; i++ )
After_Int[j][Interleaver_ST[j][i]] = Buffer[i].i;
}
void CI::Interleaver( float *input )
{
int i;
for( i=0; i<PackSize; i++ )
Buffer[i].f = input[i];
for( i=0; i<PackSize; i++ )
input[ShaffleRule[i]] = Buffer[i].f;
}
void CI::DeInterleaver( float *input )
{
int i;
for( i=0; i<PackSize; i++ )
Buffer[i].f = input[i];
for( i=0; i<PackSize; i++ )
input[i] = Buffer[ShaffleRule[i]].f;
}
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------------------
void RGCD::Init( int nt, int pathNum, int nr, int nBit, int frameLen )
{
Nt = nt;
PathNum = pathNum;
Nr = nr;
NBit = nBit;
FrameLen = frameLen;
FrameLenJ = FrameLen + PathNum - 1;
PackSize = NBit * Nt * FrameLen;
MemoryAllocation( );
switch( NBit )
{
// BPSK
case 1:
ModAng[0].Re = (float) 1; ModAng[0].Im = (float) 0;
break;
// QPSK
case 2:
ModAng[0].Re = (float) 1; ModAng[0].Im = (float) 0;
ModAng[1].Re = (float) 0; ModAng[1].Im = (float) 1;
break;
// 8PSK
case 3:
ModAng[0].Re = (float) 1; ModAng[0].Im = (float) 0;
ModAng[1].Re = (float) 0.5; ModAng[1].Im = (float) 0.866;
ModAng[2].Re = (float) -0.5; ModAng[2].Im = (float) 0.866;
break;
// 16PSK
case 4:
ModAng[0].Re = (float) 1; ModAng[0].Im = (float) 0;
ModAng[1].Re = (float) 0.7071; ModAng[1].Im = (float) 0.7071;
ModAng[2].Re = (float) 0; ModAng[2].Im = (float) 1;
ModAng[3].Re = (float) -0.7071; ModAng[3].Im = (float) 0.7071;
break;
/*
case 4:
ModAng[0].Re = (float) 1; ModAng[0].Im = (float) 0;
ModAng[1].Re = (float) 1; ModAng[1].Im = (float) 0;
ModAng[2].Re = (float) 0; ModAng[2].Im = (float) 1;
ModAng[3].Re = (float) 0; ModAng[3].Im = (float) 1;
break;
*/
}
if( NBit > 4 ) IO.ErrorMessage( "NBit > 4!!" );
}
void RGCD::MemoryAllocation( )
{
int nt, l, nr, nbit;
for( nt=0; nt<Nt; nt++ ) for( l=0; l<PathNum; l++ ) for( nr=0; nr<Nr; nr++ ) for( nbit=0; nbit<NBit; nbit++ )
{
H[nt][l][nr][nbit] = new Complex [FrameLen];
H2[nt][l][nr][nbit] = new Complex [FrameLen];
UserMean[nt][l][nr][nbit] = new Complex [FrameLen];
UserVar[nt][l][nr][nbit] = new Complex [FrameLen];
}
for( nr=0; nr<Nr; nr++ ) for( l=0; l<PathNum; l++ )
ExChipInfo[nr][l] = new float [PackSize];
for( nr=0; nr<Nr; nr++ )
{
Mean[nr] = new Complex [FrameLenJ];
Var[nr] = new Complex [FrameLenJ];
}
}
void RGCD::GetChannelInfo( char *type )
{
int nt, l, nr, nbit, len;
if( strcmp( type, "FastFade" ) == 0 )
for( nt=0; nt<Nt; nt++ ) for( l=0; l<PathNum; l++ ) for( nr=0; nr<Nr; nr++ ) for( nbit=0; nbit<NBit; nbit++ ) for( len=0; len<FrameLen; len++ )
{
H[nt][l][nr][nbit][len] = ModAng[nbit] * _H[nt][l][nr][len];
H2[nt][l][nr][nbit][len].Re = H[nt][l][nr][nbit][len].Re * H[nt][l][nr][nbit][len].Re;
H2[nt][l][nr][nbit][len].Im = H[nt][l][nr][nbit][len].Im * H[nt][l][nr][nbit][len].Im;
}
else
for( nt=0; nt<Nt; nt++ ) for( l=0; l<PathNum; l++ ) for( nr=0; nr<Nr; nr++ ) for( nbit=0; nbit<NBit; nbit++ )
{
H[nt][l][nr][nbit][0] = ModAng[nbit] * _H[nt][l][nr][0];
H2[nt][l][nr][nbit][0].Re = H[nt][l][nr][nbit][0].Re * H[nt][l][nr][nbit][0].Re;
H2[nt][l][nr][nbit][0].Im = H[nt][l][nr][nbit][0].Im * H[nt][l][nr][nbit][0].Im;
for( len=1; len<FrameLen; len++ )
{
H[nt][l][nr][nbit][len] = H[nt][l][nr][nbit][0];
H2[nt][l][nr][nbit][len] = H2[nt][l][nr][nbit][0];
}
}
}
void RGCD::Encoder( int nt, int nr, int *input, Complex *output )
{
int len, l, nbit, m, n;
for( l=0; l<PathNum; l++ ) for( n=l, m=len=0; len<FrameLen; len++, n++ ) for( nbit=0; nbit<NBit; nbit++ )
{
if( input[m] == 0 ) output[n] += H[nt][l][nr][nbit][len];
else if( input[m] == 1 ) output[n] -= H[nt][l][nr][nbit][len];
else
{
printf("ssddddddddddfffffffffffss=====%d, %d\n",m, input[m]);
IO.ImportantEorrorMessage( "Modulation error!!" );
}
m++;
}
}
void RGCD::PreDetectionSetup( char *type, float sigma2 )
{
int nt, l, nbit, len;
//for( nr=0; nr<Nr; nr++ )
for( nt=0; nt<Nt; nt++ ) for( l=0; l<PathNum; l++ ) for( nbit=0; nbit<NBit; nbit++ ) for( len=0; len<FrameLen; len++ )
{
UserMean[nt][l][0][nbit][len].Re = 0;
}
}
void RGCD::Input( float sigma2, Complex *sb, float *input )
{
float Q[LENGTH][NR2][NR2], R_ER[LENGTH][NR2];
Q_Matrix ( sigma2, &sb[0], Q, R_ER) ;
// for( i=0; i<PackSize; i++ ) input[i] = 1;
// for( nr=0; nr<Nr; nr++ )
SubInput( sigma2, &input[0], Q, R_ER );
}
void RGCD::SubInput( float sigma2, float *input, float Q[LENGTH][NR2][NR2] ,float R_ER[LENGTH][NR2])
{
int nr,nt,l, nbit,len, m, i,j, nn;
float tmp, A, tmp1, F[NR2][NR2], r_er[NR2], W[NR2],V1[NR2], HH[NR2][TOTAL_CHIP], hh[NR2];
nn=0;
for( l=0; l<PathNum; l++ ) for( len=0; len<FrameLen; len++ )
{
for( i=0; i<2*Nr; i++ )for( j=0; j<2*Nr; j++ ) F[i][j]= Q[len][i][j] ;
for( i=0; i<2*Nr; i++ )r_er[i] = R_ER[len][i];
tmp1 = QMartixFactorization_bchol(&F[0][0], 2*Nr,tmp ) ; //Q=F*F'
lowertrianglesolve( &F[0][0],&r_er[0],&W[0],2*Nr,0,2*Nr);
for( nr=0; nr<Nr; nr++ ) for( nt=0; nt<Nt; nt++ ) for( nbit=0; nbit<NBit; nbit++ )
{
m= nt*NBit +nbit;
HH[2*nr][m] = H[nt][l][nr][nbit][len].Re ; HH[2*nr+1][m] = H[nt][l][nr][nbit][len].Im ;
}
for( nt=0; nt<Nt; nt++ ) for( nbit=0; nbit<NBit; nbit++ )
{
m= nt*NBit +nbit;
for( nr=0; nr<2*Nr; nr++ )
hh[nr] = HH[nr][m];
lowertrianglesolve( &F[0][0],&hh[0],&V1[0],2*Nr,0,2*Nr);
tmp=0;tmp1=0;
for( nr=0; nr<2*Nr; nr++ )
{tmp += W[nr]*V1[nr];tmp1 += V1[nr]*V1[nr];}
A=UserMean[nt][l][0][nbit][len].Re;
A=(1- A *A ); A=(1-A*tmp1);
if (A+1.0==1.0)
IO.ImportantEorrorMessage( "Divided by O during ESE output!!!" );
tmp=2*(tmp+UserMean[nt][l][0][nbit][len].Re*tmp1)/A;
input[nn] =tmp;
nn++;
}
}
}
void RGCD::Q_Matrix (float sigma2,Complex *sb, float Q[LENGTH][NR2][NR2],float R_ER[LENGTH][NR2] )
{
int i,j ,nt, l, nr, nbit, len,m;
float matr [NR2][NR2],V[TOTAL_CHIP][TOTAL_CHIP] ,ER_tmp[NR2];
float temp[NR2][TOTAL_CHIP], HH[NR2][TOTAL_CHIP],HH_T[TOTAL_CHIP][NR2];
Complex sbsb[NR][LENGTH];
for( nr=0; nr<Nr; nr++ ) for( len=0; len<FrameLen; len++ )
sbsb[nr][len]= sb[nr*FrameLen+len];
for( l=0; l<PathNum; l++ ) for( len=0; len<FrameLen; len++ )
{
for (i=0; i<TOTAL_CHIP; i++) for (j=0; j<TOTAL_CHIP; j++ ) V[i][j] =0.0;
for( nr=0; nr<Nr; nr++ ) for( nt=0; nt<Nt; nt++ ) for( nbit=0; nbit<NBit; nbit++ )
{
m= nt*NBit +nbit;
HH[2*nr][m] = H[nt][l][nr][nbit][len].Re ; HH[2*nr+1][m] = H[nt][l][nr][nbit][len].Im ;
}
for (i=0; i<TOTAL_CHIP; i++) for (j=0; j<2*Nr; j++ ) HH_T[i][j]=HH[j][i];
m=0;
for( nt=0; nt<Nt; nt++ ) for( nbit=0; nbit<NBit; nbit++ )
{
V[m][m]=(1- UserMean[nt][l][0][nbit][len].Re *UserMean[nt][l][0][nbit][len].Re );
m++;
}
brmul(&HH[0][0],&V[0][0],2*Nr,NBit*Nt,NBit*Nt,&temp[0][0]);// Matrix product
brmul(&temp[0][0],&HH_T[0][0],2*Nr,NBit*Nt,2*Nr,&matr[0][0]);
for (i=0; i<2*Nr; i++) matr[i][i] += sigma2 ;
for (i=0; i<2*Nr; i++) for (j=0; j<2*Nr; j++) Q[len][i][j] = matr[i][j] ;
////////////////////calculate r-E(r)
for( nr=0; nr<2*Nr; nr++ ) ER_tmp[nr] = 0.0;
for( nr=0; nr<Nr; nr++ ) for( nt=0; nt<Nt; nt++ ) for( nbit=0; nbit<NBit; nbit++ )
{
ER_tmp[2*nr] += H[nt][l][nr][nbit][len].Re *UserMean[nt][l][0][nbit][len].Re;
ER_tmp[2*nr+1]+= H[nt][l][nr][nbit][len].Im *UserMean[nt][l][0][nbit][len].Re;
}
for( nr=0; nr<Nr; nr++ )
{
R_ER[len][2*nr] = sbsb[nr][len].Re -ER_tmp[2*nr];
R_ER[len][2*nr+1] = sbsb[nr][len].Im -ER_tmp[2*nr+1];
}
}
}
/*****************************************************************
function : resolve equation ax=b,
input : a----- a dim*dim lower triangular matrix
b----- a vector
dim--- dimention of the equation
begin--the elements of b from 0 to dim-1 is zero.
end ---for the future use
output : x------arguments
*****************************************************************/
int RGCD::lowertrianglesolve(float *a,float *b,float *x,int dim,int begin,int end)
{
int i,j,n;
float tb;
if(begin<0) begin = 0;
if(end>dim) end = dim;
for(i=begin;i<end;i++)
{
tb = b[i];
n = dim*i;
for(j=begin;j<i;j++)
tb -= a[n+j]*x[j];
if (a[n+i]+1.0==1.0)
{ IO.ImportantEorrorMessage( "Divided by O during the solving lower-trangle equation!!!" );
printf("fail\n"); return(-2);}
x[i]=tb/a[n+i];
}
return 1;
}
int RGCD::QMartixFactorization_bchol(float *a,int n,float det)
{ int i,j,k,u,l;
float d;
if ((a[0]+1.0==1.0)||(a[0]<0.0))
{ IO.ImportantEorrorMessage( "Divided by ZERO or <0 during the Factorization of Q!!!" );
printf("fail\n"); return(-2);}
a[0]=(float)sqrt(a[0]);
d=a[0];
for (i=1; i<=n-1; i++)
{ u=i*n; a[u]=a[u]/a[0];}
for (j=1; j<=n-1; j++)
{ l=j*n+j;
for (k=0; k<=j-1; k++)
{ u=j*n+k; a[l]=a[l]-a[u]*a[u];}
if ((a[l]+1.0==1.0)||(a[l]<0.0))
{ IO.ImportantEorrorMessage( "Divided by ZERO or <0 during the Factorization of Q_1!!!" );
printf("fail\n"); return(-2);}
a[l]=(float)sqrt(a[l]);
d=d*a[l];
for (i=j+1; i<=n-1; i++)
{ u=i*n+j;
for (k=0; k<=j-1; k++)
a[u]=a[u]-a[i*n+k]*a[j*n+k];
a[u]=a[u]/a[l];
}
}
det=d*d;
for (i=0; i<=n-2; i++)
for (j=i+1; j<=n-1; j++)
a[i*n+j]=0.0;
return(2);
}
void RGCD::brmul(float *a,float *b,int m,int n,int k,float *c)
{
int i,j,l,u;
for (i=0; i<=m-1; i++)
for (j=0; j<=k-1; j++)
{ u=i*k+j; c[u]=0.0;
for (l=0; l<=n-1; l++)
c[u]=c[u]+a[i*n+l]*b[l*k+j];
}
return;
}
void RGCD::Output( float sigma2, float *input )
{
int l, nt, nbit, len, m;
float tmp;
m=0;
// for( nr=0; nr<Nr; nr++ )
for( l=0; l<PathNum; l++ )
for( m=nt=0; nt<Nt; nt++ ) for( len=0; len<FrameLen; len++ ) for( nbit=0; nbit<NBit; nbit++ )
{
tmp= (float)exp (LClip1 ( (double) input[m] ) );
UserMean[nt][l][0][nbit][len].Re = ProbToMean( tmp );// / ExChipInfo[nr][l][m] );
m++;
}
// Average_Mean();
}
void RGCD::Average_Mean()
{
int nt, l, nr, nbit, len, s1, s2;
float Po, Ne, mini=100;
Po=0.0; Ne=0.0; s1=0; s2=0;
for( nr=0; nr<Nr; nr++ ) for( l=0; l<PathNum; l++ )
for( nt=0; nt<Nt; nt++ ) for( len=0; len<FrameLen; len++ ) for( nbit=0; nbit<NBit; nbit++ )
{
if ( UserMean[nt][l][0][nbit][len].Re >= 0 )
{
Po=Po+UserMean[nt][l][0][nbit][len].Re;
s1=s1+1;
if (UserMean[nt][l][0][nbit][len].Re <mini )
mini=UserMean[nt][l][0][nbit][len].Re;
}
else
{
Ne=Ne+UserMean[nt][l][0][nbit][len].Re;
s2=s2+1;
}
}
Po=Po/float(s1) ;Ne=Ne/float(s2) ;
for( nr=0; nr<Nr; nr++ ) for( l=0; l<PathNum; l++ )
for( nt=0; nt<Nt; nt++ ) for( len=0; len<FrameLen; len++ ) for( nbit=0; nbit<NBit; nbit++ )
{
if ( UserMean[nt][l][0][nbit][len].Re >= 0 )
{
UserMean[nt][l][0][nbit][len].Re= Po;
}
else
{
UserMean[nt][l][0][nbit][len].Re = Ne;
}
}
}
float RGCD::ProbToMean( float p )
{
return ( p - 1 ) / ( p + 1 );
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -