clear;
close all;

a = 20;
M = 10000;
X = zeros(M, 1);
for i = 1:M
    X(i) = ((a/3)^2)*randn();
end

U = zeros(M, 7);
delta = [0, 0, 0, 0, 0, 0, 0];

for R = 1:7
    quants = 2^R;
    delta(R) = 2*a/quants;
    for i = 1:M
        U(i, R) = delta(R)*fix(X(i)/delta(R));
        if X(i) >= 0
            U(i, R) = U(i, R) + delta(R)/2;
        else
            U(i, R) = U(i, R) - delta(R)/2;
        end
    end
end

SNR = zeros (7, 1);
A = X.^2;
for R = 1:7
    B = (X-U).^2;
    SNR(R) = 10*log10(sum(A)/sum(B(:,R)));
end

t = -20:0.0001:20;
f = (1/(a/3*sqrt(2*pi)))*exp(-(t.^2)/(2*((a/3)^2)));
plot(t, f);
tmp = zeros (M, 1);
epsilon = zeros (7,1);
TheorSNR = zeros (7, 1);
for R = 1:7
    for i = 1:M
        func = @(T) ((X(i)-U(i))^2)*((1/(a/3*sqrt(2*pi)))*exp(-(T.^2)/(2*((a/3)^2))));
        if X(i) >= 0
            tmp(i) = integral(func, delta(R)*fix(X(i)/delta(R)), delta(R)*fix(X(i)/delta(R))+delta(R));
        else
            tmp(i) = integral(func, delta(R)*fix(X(i)/delta(R)) - delta(R), delta(R)*fix(X(i)/delta(R)));
        end
    end
    epsilon(R) = sum(tmp);
    TheorSNR(R) = 10*log10(((a/3)^2)/epsilon(R));
end

t = 1:1:7;
figure(2);
plot(t, SNR, t, TheorSNR);