% missing sample
'---------------------------------------------------------------'

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

th = th + 0.06 * delta * randn(1, length(th));

th = th + 0.03*delta*(th .^ 2 - 1);
th = th + 0.2*delta*(th .^3 -1);
% th = th - 0.5*delta*(th .^4 -1);

% signal
K = 1;
fs = 1;
cycles = K*293+1;
% cycles = 256;
N = K*2^12;
fx = fs*cycles/N

t = linspace(0, (N-1)/fs, N);
x = sin(2*pi*fx*t) + 0.5*delta*randn(1, N);
x(198)
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:(round(fs/fx/2)-1)) * fx;
while max(fh) > 0.5
    fh = abs(fh - (fh>0.5));
end
fh(fh==0.5) = 0;
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(1);
clf;
sdb = 20*log10(s);
plot(f, sdb);
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;

figure(2);
inldnl_th(th/delta);
