## Copyright (C) 2021 Jens
## 
## 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
## <https://www.gnu.org/licenses/>.

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

## Author: Jens <jens@orion>
## Created: 2021-01-16

function iridium (name)
[x, fs] = audioread(name);
xc = (x(:,1) + i*x(:,2))';

h_rc = fir_rrc(fs, 4/fs, 0.4, 65);
xc = filter(h_rc, 1, xc);

close all;
N = length(x);

phi = angle(xc);
phi_uw = phi_unwrap(phi);
dphi = [0 diff(phi_uw)];
[lo] = losc(dphi);
xc_cfo = xc .* lo;
#xc_cfo = xc_cfo .* conj(xc_cfo);
phi_cfo_uw = phi_unwrap(angle(xc_cfo));
dphi_cfo = [0 diff(phi_cfo_uw)];
foff = dphi_cfo*fs;
found = find_preamble(dphi_cfo)

alpha_mean = 0.1;
dphi_mean = filter([0 alpha_mean], [1  -(1-alpha_mean)], dphi_cfo);

#plot(1:N, abs(xc), 1:N, atan(xc)); grid;
#plot(1:N, dphi, 1:N, phi_uw, 1:N, real(lo)); grid;

subplot(3,1,1)
#plot(1:N, dphi, 1:N, dphi_cfo); legend('dphi','dphi_{cfo}'); grid;
plot(1:N, real(xc_cfo), 1:N, imag(xc_cfo)); legend('xc_cfo_{re}','xc_cfo_{im}'); grid;

subplot(3,1,2)
plot(1:N, phi_cfo_uw, "-o"); legend('phi'); grid;

subplot(3,1,3)
plot(1:N, dphi_cfo, "-o", 1:N, dphi_mean, 'r-'); legend('dphi'); grid;

start=found(1)
ende=start + 400
step=2

figure
plot(real(xc_cfo(start:step:ende)), imag(xc_cfo(start:step:ende)), 'x'); grid;

function found = find_preamble(phi)
	thr = 0.5;
	nsmp = 120;
	last = phi(1);
	count = 0;
	found = [];
	for n=1:length(phi)
		count = count + 1;
		if abs(phi(n) - last) > thr;
			count = 0;
		endif
		if count == nsmp
			found = [found n - nsmp];
		endif
		last = phi(n);
	endfor
endfunction

function [res] = losc(dphi)
	alpha_sm = 0.995;
	dphi_sm = dphi(1);
	phi = 0;
	for n=1:length(dphi)
		res(n) = exp(-j*phi);
		phi = phi + dphi_sm;
		dphi_sm = (1-alpha_sm)*dphi_sm + alpha_sm*dphi(n);
	endfor
endfunction

function res = phi_unwrap(phi)
	offset = 0;
	last = phi(1);
	for n=1:length(phi)
		if (phi(n) - last) > pi	
			offset = offset - 2*pi;
		end
		if (phi(n) - last) < -pi	
			offset = offset + 2*pi;
		end
		res(n) = phi(n) + offset;
		last = phi(n);
	endfor
endfunction

endfunction