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
Ax = 1;                                     % input amplitude (0-peak sinusoid)
cycles = round(N * 0.1/64);                 % full-size sinewaves in simulation
fx = cycles * fs/N;                         % signal frequency
alpha = 10e-6;                              % CIN1 voltage coefficient
fx = 5e3;

[ f1, syWN1, syFS1, syt1 ] = simulate(K, N, T, 1, fx, 1e-6);
[ f2, syWN2, syFS2, syt2 ] = simulate(K, N, T, 1, fx, 10e-6);
[ f3, syWN3, syFS3, syt3 ] = simulate(K, N, T, 0.1, fx, 10e-6);
[ f4, syWN4, syFS4, syt4 ] = simulate(K, N, T, 1, 1.5*fx, 10e-6);

fixfigure(1);
plot(f1, 20*log10(syFS1+eps), 'b');  hold on;
plot(f1, 20*log10(syt1+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]  /  Int. Noise [dBV]');
legend('Output Spectrum', 'Integrated Noise');
grid on;  hold off;

fixfigure(2);
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('\alpha=10 ppm/V', '\alpha=1 ppm/V');
grid on;  hold off;

fixfigure(3);
plot(f1, 20*log10(syFS1+eps), 'b');  hold on;
plot(f3, 20*log10(syFS3+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]');
legend('A = 1V', 'A = 0.1V');
grid on;  hold off;

fixfigure(4);
plot(f1, 20*log10(syFS1+eps), 'b');  hold on;
plot(f4, 20*log10(syFS4+eps), 'r');
axis([ 0 50e3 -200 0]);
xlabel('Frequency  [Hz]');
ylabel('Output Spectrum [dBFS]');
legend(sprintf('f_{IN} = %.2gkHz', fx/1000), sprintf('f_{IN} = %.2gHz', 1.5*fx/1000));
grid on;  hold off;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f, syWN, syFS, syt] = simulate(K, N, T, 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
% 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
options = simset('Solver', 'FixedStepDiscrete');    % simulation parameter
assignin('base', 'T', T);                           % pass T to simulink

for i=1:K
    x = Ax * sin(2*pi * (fx*t + rand(1)));          % random phase
    x = x + alpha * x.^2;
    [t_, x_, Y] = sim('L22_L5_sim', max(t), options, [t', x']);
    y = Y(:, 1);
    s = abs(fft(y .* w));
    sy = sy + 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)

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