N  = 2^10;              % simulation length (output samples)
A  = 0.99;              % signal amplitude
offset = 0.00001;
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)

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

% integrate noise

figure(1);
plot(f, 20*log10(s+eps), f, 20*log10(sY+eps), 'r');
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude   [ dBWN ]');
axis([0 0.5 -80 30]);
legend('single DFT', 'average');
fixfig;
