clear all;

K  = 30;                 % number of averaged simulations
N  = 2^10;               % simulation length (output samples)
fs = 1;                  % sampling frequency 
T  = 1/fs;               % time step
fx = pi/50;              % signal frequency

A  = 0.79;               % signal amplitude
offset = pi/150;

Y = zeros(N, K);         % result vector

for i=1:K
    [t,x,y] = sim('sigma_delta_L1_sim', T*(N-1));
    Y(:, i) = y(:, 3);
end

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

figure(1);
plot(f, 10*log10(p), f, 10*log10(pint), 'r');
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude  [ dBWN ]      Power  [ dBFS ]');
axis([0 0.5 -50 30]);
legend('output spectrum [dBWN]', 'total noise [dBFS]');

% NOTE: dBFS assumes full-scale DC input, not sinusoidal input

figure(2);
plot(f, 10*log10(p+eps), f, 10*log10(pintfx+eps), 'r', f, fs./(2*f+eps), [ fs/2/16 fs/2/16 ], [ -33 16 ], 'k');
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude  [ dBWN ]      Power  [ dBFS ]');
axis([0 0.2 -50 30]);
legend('output spectrum', 'quantization noise only', 'Oversampling ratio M');
