## 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, x] = garage_ip (numSamplesPerSym)

k_n = 0.00;	% NOise
k_s = 1.0;		% Sym

baud_init = 9600;
fs = baud_init * numSamplesPerSym;

numSyms = 100;
numBursts = 3;
numSymsPauseBetweenBurst = 50;

syms = rand(numSyms, 1) > 0.5;

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

% Upconvert
x = [];

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

% Channel
[b,a] = butter(2, numSamplesPerSym/4*baud_init/fs);
y = resample (x, 131, 121);
%y = x;
y = k_s*filter(b, a, y) + k_n*randn(1, length(y));

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-5;				% Threshold forgetting factor
alpha_lead = 0.016;			% Symbol syncronizer loop filter lead 
alpha_lag = 0.000008;		% Symbol syncronizer loop filter lag 
lag_accu = 1;				% Symbol syncronizer loop filter
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 = (v_max - v_min) / 2;
	if ys < v_min
		v_min = ys;
	elseif ys > v_max
		v_max = ys;
	end
	v_min = v_min + rho_thr*(v_thr - v_min);
	v_max = v_max - rho_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
	
	% Calc error
	if (bitChg)
		err = -(count - fix(numSamplesPerSym/2));
	end
	
	% Sample bit at baud rate = fs/numSamplesPerSym
	if count == 0
		smp = [smp ys];
		smp_n = [smp_n n];
		bits = [bits bit];
	end;
	if count > 0
		count = count - 1;
	else
		count = countReload;
	end
	vcorr = alpha_lead*err + lag_accu;
	lag_accu = lag_accu + alpha_lag*err;
	baud = vcorr*fs/numSamplesPerSym;
	_err(n) = err;
	_baud(n) = baud;
	_cnt(n) = vcorr;
	_v_tr(n) = v_thr;
	_ys(n) = ys;
	n = n + 1;
	ns = ns + vcorr;
end
bits = bits';

n = 1:length(_ys);

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

endfunction