% Code for Demo. 
% One-Bit Modulo Sampling of Bandlimited Functions.
% Olga Graf, Ayush Bhandari

clc; close all; clearvars;

set(0,'defaultlinelinewidth',2)


%% Basic Function Definitions

% Following functions handles are defined that are useful later. 
% * Fractional Part
% * Rounding Function for Enforcing Amplitude Constraints
% * Pi-Bandlimited Kernel

e = exp(1); 

% Fractional Part
frac = @(x) x-floor(x);                                    

% Round Function
M = @(f,T) 2*T.*( frac(f./(2.*T) + 0.5) - 0.5 ); 

% Bandlimited Kernel
K = @(x) sinc(x);


%% Prepare Signal

% "Continuous" Time
T = 1000;
t = -10:1/T:10; t = t(:);

% Amplitudes of Random Bandlimited Signal        
Z = rand(numel(t),1) + 1i.*rand(numel(t),1);
%load Z 

% Signal
F = fft(K(t)).*Z;
f = real(ifft(F));
ampl = 9;
f = ampl*f./max(abs(f));                 % large amplitude
f = f-1;


%% Smoothing kernel

% Optimal kernel based on 2-nd order B-spline
h = 33;
B1 = (1/h)*ones(1,h);
B2 = conv(B1, B1);
B3 = conv(B2, B1);
psi = B3/max(B3);
psi = psi(1:3:end);


%% Sampling rate

% Theoretical bound on sampling rate
Tsmin = 4*pi*e*ampl*(16*(4/9)+1);

% Discrete Time
Ts = T/1;                                % sampling rate
locs = 1:T/Ts:numel(t);
ts = t(locs);
n = numel(ts);


%% Unlimited Samples 
  
% Conventional Samples   
data            = f(locs);

% Modulo Threshold
L1              = 1;

% Modulo Samples
mdata           = M(data,L1);

% Residual
resTRUE         = data-mdata;


%% Quantization 

% Quantized Conventional Samples
[u, q] = SDQ(data);

% Quantized Modulo Samples
[mu, mq] = SDQ(mdata);


%% Recovery of Residual (from Modulo Samples)

% Finite Differences of mdata
diffmdata = [zeros(1,1) diff(mdata',1)];

% Recovering the Residual
diffres = M(diffmdata,1)-diffmdata;
res = cumsum(diffres);

% Recovering data
datarec = mdata' + res;


%% Recovery of Residual (from Quantized Modulo Samples)

% Finite Differences of mq
diffmq = [zeros(1,1) diff(mq,1)];

% Finding Spikes
%
% obtain spikes (convolve fin. diff. of mq with a smoothing kernel)
spikes = conv(diffmq, psi,'same');
%
% sharpen spikes (retain only local maxima and minima)
[ymax,xmax] = findpeaks(spikes);
[ymin,xmin] = findpeaks(-spikes);
xextr = sort([xmax xmin]);
I = zeros(1,n); 
I(xextr) = 1;
sharpspikes = spikes.*I;

% Recovering the Residual
diffresSDQ = M(sharpspikes,1)-sharpspikes;
resSDQ = cumsum(diffresSDQ);

% Multibit q
qmult = mq + resSDQ;


%% Reconstruction 

% Reconstruction from Modulo Samples
frec = 0;
for k = 1:numel(ts)
    frec = frec + datarec(k).*sinc(ts-ts(k));
end
frec = frec.*mean(diff(ts));

% Reconstruction from Quantized Modulo Samples
frecSDQ = 0;
for k = 1:numel(ts)
    frecSDQ = frecSDQ + qmult(k).*sinc(ts-ts(k));
end
frecSDQ = frecSDQ.*mean(diff(ts));

% Reconstruction from Quantized Conventional Samples
rec = 0;
for k = 1:n
    rec = rec + q(k).*sinc(ts-ts(k));
end
rec = rec.*mean(diff(ts));


%% Misc

err = frec - frecSDQ;
reserr = resSDQ - resTRUE';


%% Plotting

figure('Name','')
%
subplot(2,1,1)
hold on;
axis tight
xlabel('Time')
ylabel('Amplitude (a.u.)')
plot(ts, frec, 'k'); hold on
plot(ts, frecSDQ, 'r--'); hold on
legend('Ground Truth', 'Reconstruction from One-Bit Modulo Samples')
%
subplot(2,1,2)
hold on;
axis tight
xlabel('Time')
ylabel('Amplitude (a.u.)')
plot(ts, err, 'b'); hold on
legend('Reconstruction Error')
%
FS = 12;
X = 350; Y = 220;
fig = gcf;
set(fig, 'Position', [0, 0, X, Y]);
% 
fig = gcf;
set(findall(fig,'-property','FontSize'),'FontSize',FS)
set(fig, 'Position', [0, 0, X, Y].*2.5)     
hold off;

function [u q] = SDQ(data);
    y = [data(:)];
    u = zeros(1,numel(y)+1);
    
    for k = 1:numel(y);
        q(k) = sign(u(k) + y(k) - eps);
        u(k+1) = u(k) + y(k) - q(k);
    end
end