% This matlab code compute the MSQE of the two quantization structures as in Fig 3. 
% Yuan-Pei Lin, ypl@mail.nctu.edu.tw

clear all;close all;

%generation of random variables uniform within the disc of r<=2
N0=100000;
r=2*rand(N0,1);
alpha=2*pi*rand(N0,1);
x=r.*exp(1i*alpha);

b=1:5; % the average number of quantization bits
q=2.^(b);
E1=zeros(1,size(b,2));
E1d=E1;
E2du=E1;

for i0=1:size(b,2)
    N_theta=q(i0);    
    N_phi=q(i0);     
    N_theta2=q(i0)*2;
    N_phi2=q(i0)/2;
    
%scheme1, direct quantization
theta1=(alpha+acos(r/2));
theta1=round(theta1/2/pi*N_theta)/N_theta*2*pi;
phi1=(alpha-acos(r/2));
phi1=round(phi1/2/pi*N_phi)/N_phi*2*pi;
x1=exp(1i*theta1)+exp(1i*phi1);
e1d=x1-x;

E1d(i0)=abs(e1d'*e1d)/N0;



%scheme2, direct quantization, nonuniform, noncooperated
phi2=2*acos(r/2.0000000001);
phi2=(round(phi2/pi*(N_phi2)+0.5)-0.5)/(N_phi2)*pi;
%phi2=round(phi2/pi*(N_phi2-1))/(N_phi2-1)*pi;
theta2=alpha-acos(r/2);
theta2=round(theta2/2/pi*N_theta2)/N_theta2*2*pi;
x2d=exp(1i*theta2).*(1+exp(1i*phi2));

e2dnon=x2d-x;

%scheme2, direct quantization, nonuniform, cooperated
phi2=2*acos(r/2.0000000001);
phi2=(round(phi2/pi*(N_phi2)+0.5)-0.5)/(N_phi2)*pi;
%phi2=round(phi2/pi*(N_phi2-1))/(N_phi2-1)*pi;
theta2=alpha-phi2/2;
theta2=round(theta2/2/pi*N_theta2)/N_theta2*2*pi;
x2dco=exp(1i*theta2).*(1+exp(1i*phi2));
e2dnonco=x2dco-x;


%scheme 2, direct quantization, uniform  amplitude
if N_phi2==1
    y=1;
else y=(round(r/2*(N_phi2)+.5)-0.5)/(N_phi2)*2;
end
theta2=alpha-acos(y/2);
theta2=round(theta2/2/pi*N_theta2)/N_theta2*2*pi;
x2du=exp(1i*(theta2+acos(y/2))).*y;
e2du=x2du-x;


E2du(i0)=mean(abs(e2du).^2);
E2dnon(i0)=mean(abs(e2dnon).^2);
E2dnon_co(i0)=mean(abs(e2dnonco).^2);



MSE2(i0)=pi^2/3*(4/3/N_theta2^2+1/4/N_phi2^2);
MSE2co(i0)=pi^2/3*(4/3/N_theta2^2+1/6/N_phi2^2);
MSE2u(i0)=pi^2/3*4/3/N_theta2^2+1/3/N_phi2^2;
MSE1(i0)=pi^2/3*(1/N_theta^2+1/N_phi^2);

end


%% 
figure;
semilogy(b,E1d,'b-d'); hold on,
semilogy(b, MSE1, 'b:d');
semilogy(b, E2dnon,'-rs');
semilogy(b,MSE2,':rs');
semilogy(b, E2dnon_co,'-k^');
semilogy(b,MSE2co,':k^');
semilogy(b,E2du,'k-o');
semilogy(b,MSE2u,':ko');


A=[ 'E_1, simulated                ';  
    'E_1, (4)                      ';
    'E_2, simulated                ';
    'E_2, (9)                      ';
    'E^{\prime}_2, simulated       ';
    'E^{\prime}_2, (11)            ';
    'E^{\prime\prime}_2, simulated ';
    'E^{\prime\prime}_2, (12)      ';];
legend(A);
axis([1 5 0.002 2]);
grid on
ylabel('Mean square quantization error (MSQE)');
xlabel('Average number of quantization bits b');
set(gca,'xTick',b);


