%*****************************************************************************%
%cette fonction combine les deux methodes central limite et box_muller
%on utilise 5 ROMs de 16 cases pour tabuler le LOG et un ROM de taille 256 pour le COS
%*****************************************************************************%

function gauss (Nbre_Accu,sym);

Offset_log = 0.467 
Offset_cos = 0.5
facteur_log =  1 % 511/468;
facteur_cos =  1 % 127/90;

%*****************************************************************************%
%                 Stockage des valeurs dans les ROMs
%*****************************************************************************%

for i=1:15
  rom_1(i)=floor(facteur_log*sqrt(-log((i+Offset_log)/16))*128);
  rom_2(i)=floor(facteur_log*sqrt(-log((i+Offset_log)/256))*128);
  rom_3(i)=floor(facteur_log*sqrt(-log((i+Offset_log)/4096))*128);
  rom_4(i)=floor(facteur_log*sqrt(-log((i+Offset_log)/65536))*128);
  rom_5(i)=floor(facteur_log*sqrt(-log((i+Offset_log)/1048576))*128);
end;

for i=1:256
  rom_c(i)=floor(facteur_cos*cos(pi*((i-1+Offset_cos)/512))*sqrt(2)*128);
end;

rom_1(16)=0;
rom_2(16)=0;
rom_3(16)=0;
rom_4(16)=0;
rom_5(16)=0;

%*****************************************************************************%
%                 Ecriture des valeur dans les fichiers
%*****************************************************************************%
FID=fopen('rom_cos.dat','w');
%fprintf(FID,'%5d \n', rom_1);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_c);
fclose(FID);



FID=fopen('rom_1.dat','w');
%fprintf(FID,'%5d \n', rom_1);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_1);
fclose(FID);


FID=fopen('rom_1.dat','w');
%fprintf(FID,'%5d \n', rom_1);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_1);
fclose(FID);

FID=fopen('rom_2.dat','w');
%fprintf(FID,'%5d \n', rom_2);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_2);
fclose(FID);

FID=fopen('rom_3.dat','w');
%fprintf(FID,'%5d \n', rom_3);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_3);
fclose(FID);

FID=fopen('rom_4.dat','w');
%fprintf(FID,'%5d \n', rom_4);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_4);
fclose(FID);

FID=fopen('rom_5.dat','w');
%fprintf(FID,'%5d \n', rom_5);
%fprintf(FID,'\n en Hexa : \n');
fprintf(FID,'%5X\n', rom_5);
fclose(FID);


reduction = 256

%*****************************************************************************%
%                Calcul de la distribution
%*****************************************************************************%

max = floor(1 + rom_5(1)*rom_c(1)/reduction)
epsilon = 0.000000000000001;
for i=1:max+1,
  q(i) = epsilon;
end;

for i=1:15
   for j=1:256
      x1      = floor((rom_1(i)*rom_c(j)/reduction));
      x2      = floor((rom_2(i)*rom_c(j)/reduction));
      x3      = floor((rom_3(i)*rom_c(j)/reduction));
      x4      = floor((rom_4(i)*rom_c(j)/reduction));
      x5      = floor((rom_5(i)*rom_c(j)/reduction));
      q(x1+1) = q(x1+1)+1;
      q(x2+1) = q(x2+1)+1/16;
      q(x3+1) = q(x3+1)+1/256;
      q(x4+1) = q(x4+1)+1/4096;
      q(x5+1) = q(x5+1)+1/65536;
   end;
end; 


fid=fopen('pdf_ref.dat','w');
for i = 1: max+1
   fprintf(fid,'%6.0f\n',33*q(i));
end;
status=fclose(fid);


somme =0;
box_lim=0;

% Methode de symetrisation avec biais centree sur 0.
%for i = 1 : max   
%   box(i) = 0.5*q(max+2-i);
%   box(2*max+2-i) =box(i);
%   somme  = somme + 2*box(i);
%end;
%somme=somme+q(1);
%box(max+1)=q(1);

% Methode de symetrisation avec biais centree sur -0,5.
for i = 1 : max+1   
   box(i) = 0.5*q(max+2-i);
   box(2*max+3-i) =box(i);
   somme  = somme + 2*box(i);
end;


box=box/somme;
box_lim=box;

figure(1);
subplot(211);
plot(box_lim,'g-');
grid on;

%*****************************************************************************%
%            Boucle des convolutions
%*****************************************************************************%

for i = 1 : Nbre_Accu-1
     box_lim = conv(box_lim , box);
  end;
  
subplot(212);
plot(box_lim,'g-');
grid on;
  
%*****************************************************************************%
%              Extraction de l'information relative a box_lim
%*****************************************************************************%

mC = 0;
sigma2C = 0;
sC = 0;

%*****************************************************************************%
%     Calcul des paramétres statistiques de la variable aléatoire générée
%*****************************************************************************%

for i = 1 : length(box_lim)   
   mC  = mC + i*box_lim(i);
   sigma2C  = sigma2C + i*i*box_lim(i);
end;
mC
sigma2C  = sigma2C - mC*mC;
sC = sqrt(sigma2C) 

%*****************************************************************************%
%     génération d'une variable aléatoire gaussienne ( vraie gaussienne)
%*****************************************************************************%

for i = 1: length(box_lim) 
    gauss(i) = (1/sqrt(2*pi*sigma2C))*exp(-(i-mC)*(i-mC)/(2*sigma2C));
end;

%*****************************************************************************%
%                         Représentation logarithmique
%*****************************************************************************%

logbox_lim = -log(box_lim) - 1/2*log(2*pi*sigma2C);  
loggauss   = -log(gauss) - 1/2*log(2*pi*sigma2C);


loggauss_      = abs(loggauss);     
sqrtloggauss   = sqrt(loggauss_);     

logbox_lim_    = abs(logbox_lim);
sqrtlogbox_lim = sqrt(logbox_lim_);

clear logbox_lim_; 
clear loggauss_lim_;


%*****************************************************************************%
%                        Generation courbes de reference 
%*****************************************************************************%

mCC = fix(mC); %extraction partie entiere de mC


for i=1: mCC
%  y(i)      = sqrtlogbox_lim(i+mCC-1) - sqrtloggauss(i+mCC-1);
  Epsilon(i) = (gauss(i+mCC-1) - box_lim(i+mCC-1))/gauss(i+mCC-1);
  x(i)       = (i)/sC;
end;

%*****************************************************************************%
%  Représentation graphique des erreurs : Nous faisonson une superposition
%  des courbes correspendantes aux nombre d'accumulation effectuées. Nous 
%  pouvons représenter l'erreur de réference (différence des logarithmes) 
%  ou bien l'erreur relative entre les variables aléatoires obtenues et les
%  vraies gaussienne. l'erreur est représentée entre -0.5 et 0.5, nous pouvons
%  toujours faire un Zoom pour voir cet erreur de plus prés.
%*****************************************************************************%

figure(3);
plot(x,Epsilon,'b-');
xlabel ('X/sigma')
%ylabel ('sqrt(-(Ln(P(x))+(1/2)*Ln(2*pi*sigma2)))')
ylabel ('erreur relative')
axis([0 5 -0.01 0.01])
grid on;
hold on;

%figure;
%plot(x,Err_rel)
%xlabel ('x/sigma')
%ylabel ('Erreur relative')
%axis([0 8 -0.1 0.1])
%grid on;