## 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.00);
units = struct ('x', '', 'sampleRate', '1/s', 'numSamplesPerSym', 'sps', 'loopFilterGain', '', 'with_filter', '', 'k_noise', '');

% 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;
% 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));
	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 = [];

v_min = 0;								% Threshold max
v_max = 0;								% Threshold min
hyst_thr = 0.5;							% Threshold hysteresis
rho_thr = 1e-4;							% Threshold forgetting factor
alpha_thr = 0.1;						% Min/Max update factor
k_leadLag = params.loopFilterGain;		% Overall loopfilter gain (scales k_lead and k_lag)
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
k_leak = 0;								% Lag forgetting factor
loopThreshold = 0.1;					% Processing threshold
ns = 1;									% Start sample source
n = 1;									% Start sample sink
err = 0;								% Error
while ns < N
	% Get next interpolated sample
	is = fix(ns);
	mu = ns - is;
	ys = (1-mu).*y(is) + mu.*y(is+1);

	% Calculate threshold
	v_thr = (abs(v_max) - abs(v_min)) / 2;
	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
	v_min = v_min + rho_thr*(v_thr - v_min);
	v_max = v_max + rho_thr*(v_thr - v_max);

	% Detect bit change
	bitChg = 0;
	if (bit == 1)
		if ys < (hyst_thr*v_thr)
			bit = 0;
			bitChg = 1;
		end
	else
		if ys > (hyst_thr*v_thr)
			bit = 1;
			bitChg = 1;
		end
	end
	
  isSignal = (v_thr >=loopThreshold);
	% Calc error
	err = 0;
	if (bitChg && isSignal)
		err = -(count - fix(numSamplesPerSym/2));
	end

	% Sample bit at baud rate = fs/numSamplesPerSym
	if count == 0 && isSignal
		smp = [smp ys];
		smp_n = [smp_n n];
		bits = [bits bit];
	end;
	if count > 0
		count = count - 1;
	else
		count = countReload;
	end
	vcorr = k_leadLag*k_lead*err + lag_accu;
	lag_accu = lag_accu + k_leadLag*k_lag*err;
	lag_accu = lag_accu * (1-k_leak);
	baud = fs/((1+vcorr)*numSamplesPerSym);
	_err(n) = err;
	_baud(n) = baud;
	_cnt(n) = vcorr;
	_v_tr(n) = v_thr;
	_ys(n) = ys;
	n = n + 1;
	ns = ns + (1+vcorr);
end
bits = bits';

n = 1:length(_ys);

% Plot
subplot(3, 1, 1)
plot(n, _ys, '-o', n, _v_tr, 'y', smp_n, smp, 'ro'); grid;
subplot(3, 1, 2)
plot(n, _err); grid;
subplot(3, 1, 3)
plot(n, _baud); grid;

endfunction