% specifications
fs = 1;                                     % sampling frequency
T = 1/fs;                                   % sampling period
M = 64;                                     % oversampling ratio

coeff = 'scaled';
if strcmp(coeff, 'scaled');
    % scaled filter coefficients
    a1=1;       k1=1/10;    b1=1/512;
    a2=1/2;     k2=1;       b2=1/16-1/64;
    a3=1/2;     k3=1/4;
    a4=1/4;     k4=1/4;
    a5=1/4;     k5=1/8;     
    g=2.5;
    % g = 1;
else
    % filter with rounded coefficients
    a1=1;       k1=1;       b1=1/1024;
    a2=1/2;     k2=1;       b2=1/16-1/64;
    a3=1/4;     k3=1/2;
    a4=1/8;     k4=1/4;
    a5=1/8;     k5=1/8;     g=1;
end

% simulation
K = 30;                                     % number of amplitude values
N = 2^14;                                   % simulation length
cycles = round(0.6 * N/M);                  % full-size sinewaves in simulation
fx = cycles * fs/N;                         % signal frequency
t = 0:T:(N-1)*T;                            % time base
A = logspace(-2, 0.2, K);                     % amplitude vector
options = simset('Solver', 'FixedStepDiscrete');

geff  = zeros(K);                           % result vectors
i1min = zeros(K);  i1max = zeros(K);
i2min = zeros(K);  i2max = zeros(K);
i3min = zeros(K);  i3max = zeros(K);
i4min = zeros(K);  i4max = zeros(K);
i5min = zeros(K);  i5max = zeros(K);
qmin  = zeros(K);  qmax  = zeros(K);

for i=1:K
    x = A(i) * sin(2*pi * fx*t);            % random phase sinousoid
    [t_, x_, y, i1, i2, i3, i4, i5, q] = sim('sigma_delta_L5_sim', max(t), options, [t', x']);
    geff(i) = std(y) / std(q);
    i1min(i) = min(i1);  i1max(i) = max(i1);
    i2min(i) = min(i2);  i2max(i) = max(i2);
    i3min(i) = min(i3);  i3max(i) = max(i3);
    i4min(i) = min(i4);  i4max(i) = max(i4);
    i5min(i) = min(i5);  i5max(i) = max(i5);
    qmin(i)  = min(q);   qmax(i)  = max(q);
end

Adb = 20*log10(A);

figure(1);
plot(Adb, geff, Adb, 0*Adb+0.45/k1, '-r');
xlabel('Input  [dBV]');
ylabel('Effective Quantizer Gain');
grid on;

figure(2);
color = 'bgrcmk';
for i=1:length(color)-1
    eval(sprintf('plot(Adb, i%dmin, ''%c'');', i, color(i)));  hold on;
    eval(sprintf('plot(Adb, i%dmax, ''%c'');', i, color(i)));
end
plot(Adb, qmin, 'k');
plot(Adb, qmax, 'k');
leg = [ 'i1'; 'i2'; 'i3'; 'i4'; 'i5'; 'q ' ];
for i=1:6
    y = -7-2*i;
    plot([-35 -30], [y y], color(i));
    text(-28, y, leg(i, :));
end
xlabel('Input  [dBV]');
ylabel('Loop filter peak voltages  [V]');
if strcmp(coeff, 'scaled')
    axis([min(Adb) max(Adb) -1.5 1.5]);
else 
    axis([min(Adb) max(Adb) -22 12]);
end
grid on;
hold off;
