## Copyright (C) 2018 Jens Ahrensfeld
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn {Function File} {@var{retval} =} garage (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

## Author: Jens Ahrensfeld <ahrensfeld@w2ess001vm>
## Created: 2018-12-06

function [bits] = garage_ip (varargin)

params = struct ('x', [], 'sampleRate', 48000, 'numSamplesPerSym', 16, 'loopFilterGain', 1.0, 'with_filter', 0, 'k_noise', 0.0, 'loopThreshold', 0.1);
units = struct ('x', '', 'sampleRate', '1/s', 'numSamplesPerSym', 'sps', 'loopFilterGain', '', 'with_filter', '', 'k_noise', '', 'loopThreshold', '');

% Parse parameters
names = fieldnames(params);
for k=1:nargin-1,
	for n=1:length(names),
		if strcmpi(varargin{k}, names{n})
			params.(names{n}) = varargin{k+1};
		end
	end
end
	
print_parameters(params, units)

k_n = params.k_noise;		% NOise
k_s = 1.0;		% Sym

numSamplesPerSym = params.numSamplesPerSym;
fs = params.sampleRate;
baud_init = fs / numSamplesPerSym;

numSyms = 100;
numBursts = 3;
numSymsPauseBetweenBurst = 25;

syms = rand(numSyms, 1) > 0.5;

bit = 0;
id_bit = 0;
% Line coding
syms_m = manchester(syms);

% Upconvert
if isempty(params.x)
	x = zeros(1, 2*numSymsPauseBetweenBurst*numSamplesPerSym);
	for k=1:numBursts
		for n=1:2*numSyms,
			for m=1:numSamplesPerSym,
				x = [x syms_m(n)];
			end
		end
		x = [x zeros(1, numSymsPauseBetweenBurst*numSamplesPerSym)];
	end
	x = [x zeros(1, 20*numSymsPauseBetweenBurst*numSamplesPerSym)];

	% Channel
	y = resample (x, 3200, baud_init);
	%y = x;
	if params.with_filter
		[b,a] = butter(2, numSamplesPerSym/4*baud_init/fs);
		y = k_s*filter(b, a, y) + k_n*randn(1, length(y));
	else
		y = k_s*y + k_n*randn(1, length(y));
	end
else
  x = params.x;
  y = x;
end  
N = length(y);

baud = baud_init;
bit = 0;

% The master clock counter
countReload = numSamplesPerSym-1;
count = countReload;
smp = [];
smp_n = [];
bits = [];

hyst_thr = 0.1;
v_min = 0;								% Threshold max
v_max = 0;								% Threshold min
rho_thr = 0.999;						% Threshold forgetting factor
alpha_thr = 0.01;						% Min/Max update factor
alpha_dc = 0.002;						% DC correction update factor
agc_mu = 0.01;							% AGC update factor
agc_gain_max = 10;
k_lead = 0.1;							% Symbol syncronizer loop filter lead 
k_lag = 0.0001;							% Symbol syncronizer loop filter lag 
lag_accu = 0;							% Symbol syncronizer loop filter
rho_leak = 0.9999;						% Lag forgetting factor
ns = 1;									% Start sample source
n = 1;									% Start sample sink
err = 0;								% Error

alpha_sig = 0.5/numSamplesPerSym;
v_sig = 0;

k_id = 1/numSamplesPerSym;
id_accu = 0;
dc_corr = 0;
zbit = -1;

k_lead = k_lead * params.loopFilterGain;
k_lag = k_lag * params.loopFilterGain;

isSignal = 0;

agc_target = 0.5;
agc_gain = 1.0;
while ns < N
	% Get next interpolated sample
	is = fix(ns);
	mu = ns - is;
	ys = agc_gain * ((1-mu).*y(is) + mu.*y(is+1) - dc_corr);
	
	% Calculate Runing Min/Max
	v_min = v_min * rho_thr;
	v_max = v_max * rho_thr;
	
	if ys < v_min
		v_min = (1-alpha_thr)*v_min + alpha_thr*ys;
	elseif ys > v_max
		v_max = (1-alpha_thr)*v_max + alpha_thr*ys;
	end

	% AGC
	agc_err = v_max - agc_target;
	agc_gain = max(0.0, min(agc_gain_max, agc_gain - agc_mu*agc_err));
	
	% Calculate DC offset correction
	dc_corr = (1-alpha_dc)*dc_corr + alpha_dc*(v_max + v_min);
	
	% Signal detection
	yd = abs(ys)-v_max;
	v_sig = (1-alpha_sig)*v_sig + alpha_sig*yd*yd;
	
	if isSignal == 0
		if v_sig < params.loopThreshold
			isSignal = 1;
		end
	else
		if v_sig > 1.2*params.loopThreshold
			isSignal = 0;
		end
	end	
	
	% Integrate and dump	
	id_accu = id_accu + k_id*ys;
	if count == fix(numSamplesPerSym/2)
		id_bit = id_accu >= 0.0;
		id_accu = 0;
		if isSignal
			bits = [bits id_bit];
		end
	end

	% Zero crossing detector
	bitChg = 0;
	if isSignal
		if (bit == 1)
			if ys < -hyst_thr
				bit = 0;
				bitChg = 1;
			end
		else
			if ys > +hyst_thr
				bit = 1;
				bitChg = 1;
			end
		end
	end

	% Calc error on zero crossing
	err = 0;
	if (bitChg)
		err = -(count - fix(numSamplesPerSym/2));
	end

	% Sample bit at baud rate = fs/numSamplesPerSym
	if count == 0 && isSignal
		if isSignal
			smp = [smp ys];
			smp_n = [smp_n n];
		end
	end;

	if count > 0
		count = count - 1;
	else
		count = countReload;
	end
	
	vcorr = k_lead*err + lag_accu;
	lag_accu = lag_accu * rho_leak;
	if isSignal
		if lag_accu > -0.2 && lag_accu < 0.2
			lag_accu = lag_accu + k_lag*err;
		end
	else
		lag_accu = 0;
	end
	baud = fs/((1+lag_accu)*numSamplesPerSym);
	k_id = 1/((1+lag_accu)*numSamplesPerSym);
	_err(n) = err;
	_err_baud(n) = baud - baud_init;
	_id_accu(n) = id_accu;
	_id_bit(n) = id_bit;
	_cnt(n) = vcorr;
	_v_tr(n) = agc_gain;
	_v_min(n) = v_min;
	_v_max(n) = v_max;
	_ys(n) = ys;
	n = n + 1;
	ns = max(1, ns + (1+vcorr));
end
bits = bits';

n = 1:length(_ys);

% Plot
subplot(4, 1, 1)
plot(n, _ys, '-o', smp_n, smp, 'ro'); grid;
subplot(4, 1, 2)
plot(n, _err, n, _err_baud); grid; legend('err_{smp}', 'err_{baud}')
subplot(4, 1, 3)
plot(n, _id_accu, n, _id_bit); grid; legend('accu_{id}', 'bit_{id}')
subplot(4, 1, 4)
plot(n, _v_tr, 'm', n, _v_min, n, _v_max); grid; legend('V_{sig}', 'V_{min}', 'V_{max}')

endfunction