% harmonic distortion
'---------------------------------------------------------------'

% adc with B bits and +/-1 full scale
B = 10;
delta = 2/(2^B-1);
th = -1+delta/2:delta:1-delta/2;


% signal
fs = 1;
cycles = 1501;
N = 2^14;
fx = fs*cycles/N

t = linspace(0, (N-1)/fs, N);
x = sin(2*pi*fx*t);

% nonlinear Ron
T_over_tau = 20;
Vdd_Vth = 2;
x = x .* (1 - exp(-0.5*T_over_tau ./ (1-x/Vdd_Vth)));

y = adc(x, th) * delta - 1;
s = abs(fft(y)/N*2);
s = s(1:N/2);
f = (0:length(s)-1) / N;

bx = N * fx/fs + 1;
fh = (2:4) * fx;
while max(fh) > 0.5
    fh = abs(fh - (fh>0.5));
end
bh = N * fh / fs + 1;
Asignal = s(bx);
sn = s;
sn(bx) = 1e-100;
sn(1) = 1e-100;
SFDR = 20*log10(Asignal/max(sn))
Aharm = sqrt(sum(sn(bh).^2));
sn(bh) = 1e-100;
Anoise = sqrt(sum(sn.^2));
SNR = 20*log10(Asignal/Anoise)
SNDR = 20*log10(Asignal/(Anoise+Aharm))
SDR = 20*log10(Asignal/Aharm)

figure(2);
clf;
sdb = 20*log10(s);
plot(f, sdb, 'b-');
hold on;
plot(f(bx), sdb(bx), 'ro');
plot(f(1),  sdb(1), 'y+');
plot(f(bh), sdb(bh), 'gx');
for i=1:3
    x = bh(i);
    text(f(x), sdb(x), sprintf('  H_{%d} = %.1fdBFS', i, sdb(x)));
end
text(f(1), sdb(1), sprintf('  DC = %.1fdBFS', sdb(1)));
text(f(bx), sdb(bx), sprintf('  A = %.1fdBFS', sdb(bx)));
xlabel('Frequency  [ f / f_s ]');
ylabel('Amplitude  [ dBFS ]');
title(sprintf('N = %d    SNR = %.1fdB   SDR = %.1fdB   SNDR = %.1fdB   SFDR = %.1fdB', N, SNR, SDR, SNDR, SFDR));
axis([ min(f) max(f) -120 10 ]);
fixfig;
hold off;
