hist_2_conj.cpp

source code for 2 dimensional histogram Analysis

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

using namespace std;
double PX[601][601] , PR[601][601] , X[200 * 200], X1[200 * 200], X2[200 * 200], R1[200 * 200], R2[200 * 200], R[200 * 200], DKL[50];
int M[200*200];
double GenerateWhiteGaussianNoise(double average,double sigma) ;
double quantif(double x,double Dt);
int sign(double x);
int i1 , l1 ,i2 , l2 , i3, l3 , i4 , l4, i5, l5 , i6 , l6 , i7 , l7 ; 
int m1 , n1 , k1, N , k;
double sigmaw, alpha=1.0 , p = 0.2 ;
double er, Dt, varx,L = 600;
int compt1=0 , compt2=0 ;
main()
{   
	cout << "donnez N, m et n : \n" << endl;
	cin >> N;
	cin >> m1;
	cin >> n1;

	//Initialisation 


for (k=0;k<40;k++)
{
	DKL[k]=0.0;
}

for (k1=0;k1<40;k1++)
{
		cout << k1 << endl;
		for(i1=0;i1<=L;i1++)
		{
			for(l1=0;l1<=L;l1++)
			{   
				PX[i1][l1]=0.0 ;
				PR[i1][l1]=0.0 ;
			}
		}
	for(i2=0;i2<N;i2++)
	{

		//Generation du signal gaussien
		for(i3=0;i3 < m1*n1-1;i3++)
		{

			double interm;

			{   sigmaw=20*10^(-k1/10);
				X[i3] = GenerateWhiteGaussianNoise(0.0, 20) ;
				interm = GenerateWhiteGaussianNoise(0.0, 1.0);
				M[i3]=(sign(interm)+1)/2;
				Dt=sqrt(12*sigmaw);
				er=quantif(X[i3]-M[i3]*Dt/2,Dt)-(X[i3]-M[i3]*Dt/2);
				R[i3]=X[i3]+alpha*er;

					if (i3 %2==0)
					{

					X1[compt1] = X[i3];
					R1[compt1] = R[i3];
					compt1=compt1++;
cout << compt1 <<endl;

					}
					else 
					{

					X2[compt2] = X[i3];
					R2[compt2] = R[i3];
					compt2=compt2++;
cout << compt2 <<endl;
					} 

            
			}
		}


	//Calcul de l'histograme
	for(i4=0;i4<L;i4++)
	{  
		for(l4=0;l4<L;l4++)
		{ 
			for(i5=0;i5<m1*n1-1;i5++)
			{  

				    if (((i4-1-L/2)*p<=X1[i5])&&(X1[i5]<(i4-L/2)*p)&&((l4-1-L/2)*p<=X2[i5])&&(X2[i5]<(l4-L/2)*p))
					{
						PX[i4][l4]	=	PX[i4][l4] + 2.0/(m1*n1*p);
					}

				    if (((i4-1-L/2)*p<=R1[i5])&&(R1[i5]<(i4-L/2)*p)&&((l4-1-L/2)*p<=R2[i5])&&(R2[i5]<(l4-L/2)*p))
					{
						PR[i4][l4]	=	PR[i4][l4] + 2.0/(m1*n1*p);
//						PR[i4+601*l4]	=	PR[i4+601*l4] + 2.0/(m1);
					}

			}
		}
	}
	//Calcul de la KLD_2D	
			for(l6=0;l6<=L;l6++)
			{
				if ((PX[i4][l6]!=0.0)&&(PR[i4][l6]!=0.0))
				{
					DKL[k1]=DKL[k1]+PX[i4][l6]*log(PX[i4][l6]/PR[i4][l6])/log(2.0);
				}
			}

	//Affichage des r閟ultats 
//	for(i7=0;i7<L;i7++)
//	{  
//		for(l7=0;l7<L;l7++)
//		{
//                    if( PX[i7][l7]!=1.0e-21 )
//					{
//					cout << "PR[" <<i7<< "]["<<l7<<"]=" << PR[i7][l7] << "\n PX["<<i7<<"]["<<l7<<"]=" << PX[i7][l7]<<endl;
//					}
//		}
//	}
	}
	cout << DKL[k1] << endl;
}
return 0;
}


//*******Fonction qui genere l'echantillon gaussien*******
double GenerateWhiteGaussianNoise(double average,double sigma) 
{
  double M_PI=3.14;
  double x,y;
  double g,h;
   do {
     x = (double)rand() / RAND_MAX;
     y = (double)rand() / RAND_MAX;
  } while(x == 0.0 || y == 0.0 || x == 1.0 || y == 1.0);

  g = sqrt(-log(x));
  h = sqrt(2.0) * cos(2.0 * M_PI * y);
   return (average + sigma * g * h);
}

//****************Fonction de quantification**********************
	double quantif(double x,double Dt)
	{double Q;
     Q=Dt*floor(x/Dt+0.5);
	 return Q;
	}
//**************fonction calcul sign***************************
int sign(double x)
{int y;
if (x==0)
	y=0;
else 
	{if (x>0) 
	     y=1;
	 else
		 y=-1;
	}
	return y;
}