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

% 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
g = 3;                                      % feedback gain
A  = 0;                                     % 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
x = 0.006 * ones(1, N);                     % input (offset)
w = hodiewindow(N);                         % window
syw = zeros(N/2, 1);                        % averaged windowed output spectrum
options = simset('Solver', 'FixedStepDiscrete');

Afs2_ = [ 0, 1e-6, 1e-3 ];
syw = zeros(N/2, length(Afs2_));
iyw = syw;
for j=1:length(Afs2_)
    Afs2 = Afs2_(j)
    for i=1:K
        x(1) = rand(1) - 0.5;                   % random first term to enable averaging
        [t_, x_, y, i1, i2, i3, i4, i5, q] = sim('L21_L5_Vref_sim', max(t), options, [t', x']);
        s = abs(fft(y.*w));
        syw(:,j) = syw(:,j) + s(1:N/2);
    end

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

end

% display results

ax = [ 0          50e3 -150 50; ...
       0          fs/2 -150 50; ...
       0.49*fs    fs/2 -150 50 ];

st = [ 'b'; 'c'; 'g'; 'm'; 'k' ];
        
for i=1:size(ax, 1)
    fixfigure(i);
    leg = cell(length(Afs2_), 1);
    for j=1:length(Afs2_)
        plot(f, 20*log10(syw(:,j)+eps), st(j), 'linewidth', 1.8);  
        hold on;
        leg(j) = cellstr(sprintf('%gV', Afs2_(j)));
    end
    xlabel('Frequency  [Hz]');
    ylabel('Output Spectrum [dBWN]');
    legend(leg, 0);
    axis(ax(i,:));
    grid on;  hold off;
end
