% specifications
fs = 3e6;                                   % 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;

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

% simulation
K = 30;                                     % number of simulation averages
N = 2^15;                                   % simulation length
A  = 0.3;                                   % signal amplitude
cycles = round(0.1 * N/M);                  % full-size sinewaves in simulation
fx = cycles * fs/N;                         % signal frequency

f = linspace(0, fs/2, N/2);
t = 0:T:(N-1)*T;                            % time base
w = hodiewindow(N);                         % window
syw = zeros(N/2, 1);                        % averaged windowed output spectrum
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('L20_L5_sim', max(t), options, [t', x']);
    s = abs(fft(y.*w));
    syw = syw + s(1:N/2);
end

syw = syw/sqrt(N)/std(w)/K;
% syw(1:5) = 0;                               % chop off DC
iyw = sqrt(cumsum(syw.^2)/N);

% display results

fixfigure(1);
plot(f, 20*log10(syw+eps), 'b', 'linewidth', 1.5);  hold on;
plot(f, 20*log10(iyw+eps), 'r');
axis([0 50e3 -150 50]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBV]');
ylabel('Output Spectrum [dBWN]');
legend('Output Spectrum', sprintf('Integrated Noise (%d averages)', K), 4);
grid on;  hold off;

fixfigure(2);
plot(f, 20*log10(syw+eps), 'b', 'linewidth', 1.5);  hold on;
plot(f, 20*log10(iyw+eps), 'r');
axis([0 fs/2 -150 50]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBV]');
ylabel('Output Spectrum [dBWN]');
legend('Output Spectrum', sprintf('Integrated Noise (%d averages)', K), 4);
grid on;  hold off;

fixfigure(3);
plot(f, 20*log10(syw+eps), 'b', 'linewidth', 1.5);  hold on;
plot(f, 20*log10(iyw+eps), 'r');
axis([0.49*fs 0.5*fs -150 50]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBV]');
legend('Output Spectrum', sprintf('Integrated Noise (%d averages)', K), 4);
grid on;  hold off;

