From 7293ec803f51ae1baf30f8c86444f546a0e7b0e7 Mon Sep 17 00:00:00 2001
From: Jens Ahrensfeld <jens@jayfield.org>
Date: Mon, 7 Jan 2019 14:54:30 +0000
Subject: [PATCH] - improvement: use Integrate and Dump filter

git-svn-id: http://moon:8086/svn/matlab/trunk@113 801c6759-fa7c-4059-a304-17956f83a07c
---
 garage_ip.m | 70 +++++++++++++++++++++++++++++++++--------------------
 1 file changed, 44 insertions(+), 26 deletions(-)

diff --git a/garage_ip.m b/garage_ip.m
index c006b40..47cbdab 100755
--- a/garage_ip.m
+++ b/garage_ip.m
@@ -24,7 +24,7 @@
 
 function [bits] = garage_ip (varargin)
 
-params = struct ('x', [], 'sampleRate', 48000, 'numSamplesPerSym', 16, 'loopFilterGain', 1.0, 'with_filter', 0, 'k_noise', 0.00);
+params = struct ('x', [], 'sampleRate', 48000, 'numSamplesPerSym', 16, 'loopFilterGain', 1.0, 'with_filter', 0, 'k_noise', 0.0);
 units = struct ('x', '', 'sampleRate', '1/s', 'numSamplesPerSym', 'sps', 'loopFilterGain', '', 'with_filter', '', 'k_noise', '');
 
 % Parse parameters
@@ -53,6 +53,7 @@ numSymsPauseBetweenBurst = 25;
 syms = rand(numSyms, 1) > 0.5;
 
 bit = 0;
+id_bit = 0;
 % Line coding
 syms_m = manchester(syms);
 
@@ -75,6 +76,8 @@ if isempty(params.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;
@@ -92,65 +95,77 @@ smp = [];
 smp_n = [];
 bits = [];
 
+hyst_thr = 0.1;
 v_min = 0;								% Threshold max
 v_max = 0;								% Threshold min
-hyst_thr = 0.5;							% Threshold hysteresis
-rho_thr = 1e-3;							% Threshold forgetting factor
-alpha_thr = 0.1;						% Min/Max update factor
+rho_thr = 0.999;						% Threshold forgetting factor
+alpha_thr = 0.002;						% 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.05;					% Processing threshold
+k_leak = 0.0;							% Lag forgetting factor
+loopThreshold = 0.02;					% Processing threshold
 ns = 1;									% Start sample source
 n = 1;									% Start sample sink
 err = 0;								% Error
 
 alpha_sig = 0.5/numSamplesPerSym;
-alpha_sig_m = alpha_sig/5;
 v_sig = 0;
-v_sig_m = 0.5;
 
+k_id = 1/numSamplesPerSym;
+id_accu = 0;
+dc_corr = 0;
+zbit = -1;
 while ns < N
 	% Get next interpolated sample
 	is = fix(ns);
 	mu = ns - is;
-	ys = (1-mu).*y(is) + mu.*y(is+1);
+	ys = (1-mu).*y(is) + mu.*y(is+1) - dc_corr;
 
-	% Calculate threshold
-	v_thr = (abs(v_max) - abs(v_min)) / 2;
+	% 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;
+		v_min = ys;
 	elseif ys > v_max
-		v_max = (1-alpha_thr)*v_max + alpha_thr*ys;
+		v_max = ys;
 	end
-	v_min = v_min + rho_thr*(v_thr - v_min);
-	v_max = v_max + rho_thr*(v_thr - v_max);
 
+	% Calculate DC offset correction
+	dc_corr = dc_corr + alpha_thr*(v_max + v_min);
+	
 	% Signal detection
-	y_nodc = ys-v_thr;
-	v_sig = (1-alpha_sig)*v_sig + alpha_sig*y_nodc*y_nodc;
-	v_sig_m = (1-alpha_sig_m)*v_sig_m + alpha_sig_m*ys;
+	v_sig = (1-alpha_sig)*v_sig + alpha_sig*ys*ys;
 	isSignal = v_sig >= loopThreshold;
 	
-	% Detect bit change
+	% 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*v_thr)
+			if ys < -hyst_thr
 				bit = 0;
 				bitChg = 1;
 			end
 		else
-			if ys > (hyst_thr*v_thr)
+			if ys > +hyst_thr
 				bit = 1;
 				bitChg = 1;
 			end
 		end
 	end
 
-	% Calc error
+	% Calc error on zero crossing
 	err = 0;
 	if (bitChg)
 		err = -(count - fix(numSamplesPerSym/2));
@@ -158,12 +173,12 @@ while ns < N
 
 	% Sample bit at baud rate = fs/numSamplesPerSym
 	if count == 0 && isSignal
-		bits = [bits bit];
 		if isSignal
 			smp = [smp ys];
 			smp_n = [smp_n n];
 		end
 	end;
+
 	if count > 0
 		count = count - 1;
 	else
@@ -173,8 +188,11 @@ while ns < N
 	lag_accu = lag_accu + k_leadLag*k_lag*err;
 	lag_accu = lag_accu * (1-k_leak);
 	baud = fs/((1+lag_accu)*numSamplesPerSym);
+	k_id = 1/((1+lag_accu)*numSamplesPerSym);
 	_err(n) = err;
-	_baud(n) = baud;
+	_err_baud(n) = baud - baud_init;
+	_id_accu(n) = id_accu;
+	_id_bit(n) = id_bit;
 	_cnt(n) = vcorr;
 	_v_tr(n) = v_sig;
 	_ys(n) = ys;
@@ -189,8 +207,8 @@ n = 1:length(_ys);
 subplot(3, 1, 1)
 plot(n, _ys, '-o', n, _v_tr, 'm', smp_n, smp, 'ro'); grid;
 subplot(3, 1, 2)
-plot(n, _err); grid;
+plot(n, _err, n, _err_baud); grid; legend('err_{smp}', 'err_{baud}')
 subplot(3, 1, 3)
-plot(n, _baud); grid;
+plot(n, _id_accu, n, _id_bit); grid; legend('accu_{id}', 'bit_{id}')
 
 endfunction