%**********************************************************
% Emmanuel Boutillon               
% Matlab file                            
% LESTER, Universite de Bretagne Sud , July 21th 2000
%**********************************************************
function wgng(Nbre_Accu);

K        = 5   %5
q_f      = 4   %4
q_g      = 8   %8
m_f      = 7   %7   % 3+m_f bits to represent f(x)>0 3.m_f      
m_g      = 6   %6  % 1+m_f bits to represnet g(x)   1.m_g 
d_f  = 0.36 %0.36   
d_g  = 0.5  %0.5
b = 6

%*** Creation of the ROMs fr(x) ***************%
for r=1:K
  for s=1:pow2(q_f)-1  
    rom_f(s,r)=floor(sqrt(-log((s+d_f)*pow2(-r*q_f)))*pow2(m_f));
  end;
  rom_f(pow2(q_f), r) = 0;
end;

%*** Creation of the ROM g(x) ***************%
for u=1:pow2(q_g)
  rom_g(u)=floor(sqrt(2)*cos(pi*((u-1+d_g)*pow2(-q_g-1)))*pow2(m_g));
end;

scaling = pow2(m_f + m_g - b);  
% To have exactly b bits after the dot

%*** Initialisation of the distribution HBM
max = floor(1 + rom_f(1,K)*rom_g(1)/scaling)
HBM = zeros(1,max+1);

% Construction of half the distribution
for s=1:pow2(q_f)-1
  for r=1:K
    for u=1:pow2(q_g)
   n=floor((rom_f(s,r)*rom_g(u)/scaling));
   HBM(n+1) = HBM(n+1) + pow2(-(r*q_f + q_g));
   end; 
 end; 
end; 

% Symetrisation of the distribution
for i = 1 : max   
   BM(i) = 0.5*HBM(max+2-i);
   BM(2*max+2-i) = BM(i);
end;
BM(max+1)=HBM(1);  % probability of n=0

%*** Convolution loop
BMi = BM;
for i = 1 : Nbre_Accu-1
     BMi = conv(BMi , BM);
end;    

box_lim = BMi;

%%
%*********************************************************
%  means and standard deviation
%*********************************************************
%%
mC = 0;
sigma2C = 0;
sC = 0;

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


%*****************************************************************************%
%   Generation of N(mC, SC)
%*****************************************************************************%
for i = 1: length(box_lim) 
   gauss(i) = (1/sqrt(2*pi*sigma2C))*exp(-(i-mC)*(i-mC)/(2*sigma2C));
   xx(i) = (i-mC)/sC;
end;
figure(2);
plot(xx,box_lim,'b',xx,gauss,'r');

%*****************************************************************************%
%  Function epsilon 
%*****************************************************************************%
mCC = fix(mC); %integer value of mC

for i=1: mCC
  x(i)       = (i)/sC;
  epsilon(i) = (gauss(i+mCC-1) - box_lim(i+mCC-1))/gauss(i+mCC-1);
end;

%*****************************************************************************%
%  Plot curves
%*****************************************************************************%
figure(1)
plot(x,epsilon);
xlabel ('x')
ylabel ('Epsilon(x)')
axis([0 5 -0.1 0.1])
grid on;
hold on;

