function L22_L5
% 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
[ f, syWN, syFS, syt ] = simulate(K, N, T, Ax, fx);

% output
plotFS(1, f, syFS, syt, 0, 50e3);           % dBFS
plotFS(2, f, syFS, syt, fs/2-50e3, fs/2);
plotFS(3, f, syFS, syt, 0, fs/2);

plotWN(4, f, syWN, syt, 0, 50e3);           % dBWN
plotWN(5, f, syWN, syt, fs/2-50e3, fs/2);
plotWN(6, f, syWN, syt, 0, fs/2);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotFS(fig, f, sy, syt, fmin, fmax)
% plot output spectrum
% fig figure number
% f frequency
% sy spectral density
% syt integrated noise
% fmin, fmax x-axis range

fixfigure(fig);
plot(f, 20*log10(sy+eps), 'b', 'linewidth', 1.5);  hold on;
if length(syt) == length(f)
    plot(f, 20*log10(syt+eps), 'r');
    ylabel('Output Spectrum [dBFS]  /  Int. Noise [dBV]');
    legend('Output Spectrum', 'Integrated Noise', 4);
else
    ylabel('Output Spectrum [dBFS]');
end
axis([ fmin fmax -200 0]);
xlabel('Frequency  [Hz]');
grid on;  hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotWN(fig, f, sy, syt, fmin, fmax)
% plot output spectrum
% fig figure number
% f frequency
% sy spectral density
% syt integrated noise
% fmin, fmax x-axis range

fixfigure(fig);
plot(f, 20*log10(sy+eps), 'b', 'linewidth', 1.5);  hold on;
if length(syt) == length(f)
    plot(f, 20*log10(syt+eps), 'r');
    ylabel('Output Spectrum [dBWN]  /  Int. Noise [dBV]');
    legend('Output Spectrum', 'Integrated Noise', 4);
else
    ylabel('Output Spectrum [dBFS]');
end
axis([ fmin fmax -200 50]);
xlabel('Frequency  [Hz]');
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