N  = 2^10;              % simulation length (output samples)
A  = 0.99;              % signal amplitude
offset = 0.001;
fs = 1;                 % sampling frequency 
T  = 1/fs;              % time step
fx = 0.1 + 0.000001;    % signal frequency
sigma = 1;

[t,x,y] = sim('L19_level1_dither_sim', T*(N-1));

X = y(:,1);
Q = y(:,2);
Y = y(:,3);
w = nuttallwin(length(Y));

s = abs(fft(Y .* w)) / sqrt(N);
s = s(1:N/2);
f = (0:length(s)-1) / N;
G = std(Y)/std(Q)

figure(1);
% average K spectra
sY = zeros(N/2, 1);
K = 30;
for k=1:K
    offset = 0.005*rand(1);
    [t,x,y] = sim('L19_level1_dither_sim', T*(N-1));
    Y = y(:,3);
    sy = abs(fft(Y.*w/sqrt(N)));
    sY = sY + sy(1:N/2);
end
sY = sY / K / std(w);
sY(1) = eps;

int = cumsum(sY.*sY) / N;

% integrate noise

figure(1);
plot(f, 20*log10(sY+eps), 'r', f, 10*log10(int+eps));
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude  [ dBWN ]   /    Integrated Noise  [dBV]');
legend('Output Spectrum', 'Integrated Noise');
axis([0 0.5 -80 30]);
fixfig;
