%函數quantize.m:高斯概率密度函數的非均勻量化%
function [r,d]=quantize(L)                %L:重構層數         
r=linspace(0,0,L); d=linspace(0,0,L+1);
d(1)=0; d(L+1)=500;                       %設定d0=0，dL =500(代表無窮大)
d_inc=0.005; r_inc=0.003; error_bound1=0.002; error_bound2=0.003;  
r_estimate=500;tolerance=[1e-4 1e-4];                        
r_init=0.01;                              %r(1)值   
check=0;
while check==0
  r(1)=r_init;                
  for k=1:L-1
    d(k+1)=d(k);                          %猜測的d(k+1)以d(k)開始運算          
    %比對所得的r_estimate和r(k)之差，此循環在r(k)給定后找出d(k+1)% 
    while abs(r_estimate-r(k))>error_bound1       
      d(k+1)=d(k+1)+d_inc;                %驗算不符合略微增加d_inc
      Int_N=quad('gau1',d(k),d(k+1),tolerance); 
      Int_D=quad('gau',d(k),d(k+1),tolerance); 
      r_estimate=Int_N/Int_D;             
    end
    r(k+1)=2*d(k+1)-r(k);                 %從d(k+1)與r(k)求r(k+1)
  end                                      
  %從公式算出的rL_estimate和r(k+1)做比較% 
  IntL_N=quad('gau1',d(L),d(L+1),tolerance);
  IntL_D=quad('gau',d(L),d(L+1),tolerance); 
  rL_estimate=IntL_N/IntL_D;
  diff=rL_estimate-r(L);
  relative_error=abs(diff)/rL_estimate;                                       
        
   if relative_error<error_bound2
    check=1;
  else
    check=0;
    r_init=r_init+r_inc;                  %不符合條件則重新選取r(1)
  end        
end