Nfft = 4096;
fs = 245.76e6;
f1 = 13.5e6;
f2 = 15.3e6;

t = 0:1/fs:(Nfft-1)/fs;
vin = sin(2*pi*f1*t)/2+sin(2*pi*f2*t)/2;
%9.2bit thermal noise
vin = vin + randn(1, Nfft) * 1 / 2^9.2 / 2;

%Ideal 10bit quantization
out = round(vin*2^10/2)+2^9;

f = linspace(0, 0.5, Nfft/2+1)*fs;
OUT = fft(out, Nfft);
figure(1);
plot(f, 20*log10(abs(OUT(1:Nfft/2+1))));

idx = round(f1/fs*Nfft)+1;
idx2 = round(f2/fs*Nfft)+1;

sig_power = abs(OUT(idx)).^2 + abs(OUT(idx2)).^2;
noise_power = sum(abs(OUT(2:Nfft/2)).^2)-sig_power;
SNR = 10*log10(sig_power/noise_power)+3.01;
display(sprintf('SNR: %.2fdB, ENOB: %.2f bits', SNR, (SNR-1.76)/6.02));

%kill two samples
out2 = out;
out2(7:8:end)=0;
out2(8:8:end)=0;

OUT2 = fft(out2, Nfft);

figure(2);
plot(f, 20*log10(abs(OUT2(1:Nfft/2+1))));

%Find the "mixing" function
mix = out2~=0;
MIX = fft(mix, Nfft);
MIXR = fliplr(MIX);

%Create a sparse matrix whose entries correspond to the reversed and
%shifted FFT of the mixing function, as used in the periodic convolution
num_impulses = sum(MIX~=0);
r = filter(ones(1, num_impulses), 1, upsample(1:Nfft, num_impulses));
c = zeros(Nfft, num_impulses);
v = zeros(Nfft, num_impulses);

row = zeros(1,Nfft);
%Go row by row to deal with the circular shift; this is slow and could
%potentially be optimized
for m=1:Nfft
    row = circshift(MIXR, m, 2);
    c(m, :) = find(row ~= 0);
    v(m, :) = row(c(m, :));
end
c = reshape(c', 1, Nfft*num_impulses);
v = reshape(v', 1, Nfft*num_impulses);

Mk = sparse(r, c, v);

OUT_GUESS = abs(Mk)\(Nfft*abs(OUT2)');

figure(3);
plot(f, 20*log10(abs(OUT_GUESS(1:Nfft/2+1))));

idx = round(f1/fs*Nfft)+1;
idx2 = round(f2/fs*Nfft)+1;

display('Guessed');

sig_power = abs(OUT_GUESS(idx)).^2 + abs(OUT_GUESS(idx2)).^2;
noise_power = sum(abs(OUT_GUESS(2:Nfft/2)).^2)-sig_power;
SNR = 10*log10(sig_power/noise_power)+3.01;
display(sprintf('SNR: %.2fdB, ENOB: %.2f bits', SNR, (SNR-1.76)/6.02));