function L22_L5_nonlinear_CIN1
% specifications
fs = 3e6;                                   % sampling frequency
T = 1/fs;                                   % sampling period

% simulation
K = 30;                                      % number of simulation averages
N = 2^15;                                   % simulation length
fx = 5e3;;                                  % signal frequency

alpha1 = 10e-6;
alpha2 = 1000e-6;

[ f1, syWN1, syFS1, syt1, si1FS1 ] = simulate(K, N, T, 0, 1, fx, alpha1);
[ f2, syWN2, syFS2, syt2, si1FS2 ] = simulate(K, N, T, 0, 1, fx, alpha2);
[ f3, syWN3, syFS3, syt3, si1FS3 ] = simulate(K, N, T, 0.005, 0, fx, alpha2);
[ f4, syWN4, syFS4, syt4, si1FS4 ] = simulate(K, N, T, 0.008, 0, fx, alpha2);
[ f5, syWN5, syFS5, syt5, si1FS5 ] = simulate(K, N, T, 0.005, 0, fx, alpha1);

fixfigure(1);
plot(f1, 20*log10(syFS1+eps), 'b');  hold on;
plot(f2, 20*log10(syFS2+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]');
legend(sprintf('\\alpha_{FB1} = %g ppm/V', alpha1*1e6), sprintf('\\alpha_{FB1} = %g ppm/V', alpha2*1e6))
grid on;  hold off;

fixfigure(2);
plot(f1, 20*log10(syFS1+eps), 'b');  hold on;
plot(f1, 20*log10(si1FS1+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]');
legend('Modulator Output', 'Output of 1st Integrator');
grid on;  hold off;

fixfigure(3);
plot(f3, 20*log10(syFS3+eps), 'b');  hold on;
plot(f4, 20*log10(syFS4+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]');
legend('5mV DC', '8mV DC');
grid on;  hold off;

fixfigure(4);
subplot(2, 1, 1);
plot(f3, 20*log10(syFS3+eps), 'b');
ylabel('Output Spectrum [dBFS]');
axis([ 0 50e3 -200 0]);
grid on;
subplot(2, 1, 2);
plot(f3, 20*log10(syFS3+eps), 'r');
ylabel('Output Spectrum [dBFS]');
axis([ fs/2-50e3 fs/2 -200 0]);
xlabel('Frequency  [Hz]');
grid on;

fixfigure(5);
subplot(2, 1, 1);
plot(f5, 20*log10(syFS5+eps), 'b');
ylabel('Output [dBFS]', 'FontSize', 14);
axis([ fs/2-50e3 fs/2 -200 0]);
ha = gca;
set(ha, 'FontSize', 14);
set(ha, 'FontName', 'Arial');
set(ha, 'LineWidth', 2);
grid on;
subplot(2, 1, 2);
plot(f5, 20*log10(syFS5+eps), 'r');
ylabel('Output [dBFS]');
xlabel('Frequency  [Hz]');
axis([ 0 50e3 -200 0]);
ha = gca;
set(ha, 'FontSize', 14);
set(ha, 'FontName', 'Arial');
set(ha, 'LineWidth', 2);
grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f, syWN, syFS, syt, si1FS] = simulate(K, N, T, Adc, Ax, fx, alpha)
% simulate 5th order modulator
% K number of runs to average
% N number of samples (over which DFT is run)
% T sampling time
% Adc offset
% Ax amplitude of (sinusoidal) input
% fx frequency of (sinusoidal) input

% scaled filter coefficients
a1=1;       k1=1/10;    b1=1/512;           g=3;
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;

w = hodiewindow(N);                                 % window
t = 0:T:(N-1)*T;                                    % time vector
f = linspace(0, 0.5/T, N/2);                        % frequency vector
sy = zeros(N/2, 1);                                 % averaged windowed output spectrum
si1= sy;
options = simset('Solver', 'FixedStepDiscrete');    % simulation parameter
assignin('base', 'T', T);                           % pass T to simulink
assignin('base', 'alpha', alpha);

for i=1:K
    x = Ax * sin(2*pi * (fx*t + rand(1)));          % random phase
    x = x + Adc * ones(1, N);
    x(1) = rand(1) - 0.5;                           % randomize dc simulations
    % x = x + alpha * x.^2;                           % Cin1 nonlinearity
    [t_, x_, Y] = sim('L22_L5_nonlinear_CFB_sim', max(t), options, [t', x']);
    y = Y(:, 1);
    s = abs(fft(y .* w));
    sy = sy + s(1:N/2);
    i1 = Y(:, 2);                                   % output from 1st integrator
    s = abs(fft(i1 .* w));
    si1 = si1 + s(1:N/2);
end

sy = sy / K;                                        % scale:
syWN = sy / sqrt(N);                                %    relative to white noise
syFS = sy / N * g * 2;                              %    relative to full scale (1V 0-peak)
si1FS= si1 / N * g * 2;

syt = syWN;                                         % total noise (integrates to 1 iff signal is not excluded)
syt(1:7) = 0;                                       % exclude DC bins
nx = round(N * fx * T + 1);                         % fundamental center bin
syt(nx-7:nx+7) = 0;                                 % exclude fundamental
syt = sqrt(cumsum(syt .^2)/length(syt));            % noise integral
