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

% 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;

k = [ k1 k2 k3 k4 k5 ];
a = [ a1 a2 a3 a4 a5 ];
b = [ b1 b2 ];

% simulation
K = 20;                                     % number of simulation averages
N = 2^15;                                   % simulation length
A  = 0.056;                                 % signal amplitude
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
Y = zeros(N, K);                            % result vector
options = simset('Solver', 'FixedStepDiscrete');

for i=1:K
    x = A * sin(2*pi * (fx*t + rand(1)));   % random phase sinousoid
    % x = 0.002 * ones(1, N);               % dc
    % x(1) = rand(1) - 0.5;                 % random first term to enable averaging
    % x = 9e-6 * 2 * (rand(1, N)-0.5) + x;  % dither at input
    [t_, x_, y, i1, i2, i3, i4, i5, q] = sim('sigma_delta_L5_sim', max(t), options, [t', x']);
    Y(:, i) = y;
end

[f, p, pint, pintfx] = avg_spectrum(Y, fx/fs);

figure(1);
plot(f, 10*log10(p+eps), f, 10*log10(pint+eps));
xlabel('Frequency  [f/fs]');
ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBFS]');
legend('Output Spectrum', sprintf('Integrated Noise (%d averages)', K), 4);
grid on;
axis([0 0.5*fs -165 45]);

figure(2);
plot(f*M, 10*log10(p+eps), f*M, 10*log10(pint+eps));
xlabel('Frequency  [f/fN]');
ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBFS]');
legend('Output Spectrum', sprintf('Integrated Noise (%d averages)', K), 2);
grid on;
axis([0 fs -165 45]);
