## Copyright (C) 2018 -*- octave -*-

## usage: [x,xx,fs,last_gpsfix]=proc_kiwi_iq_wav(wav_filename, max_last_gpsfix)
##
## input: wav_filename - filename of .wav file generated by kiwirecorder.py with options '-m iq --kiwi-wav'
##        (optional) 2nd argument - maximal allowed age of the GNSS position solution (default is max_last_gpsfix=255)
## output: x - array of structs containing the information in the WAV chunks, i.e. IQ samples and GNSS timestamps
##         xx - array of structs containting IQ samples and for each IQ samples a GNSS timestamp
##         fs - sampling rate
##         last_gpsfix - oldest GNSS position solution
##
## note: last_gps_fix = 255 means there was no GNSS position solution at all on the KiwiSDR

function [x,xx,fs,last_gpsfix]=proc_kiwi_iq_wav(fn, max_last_gpsfix)
  if nargin == 0
    print_usage;
  end
  if nargin == 1
    max_last_gpsfix = 255;
  end

  [x,fs]      = read_kiwi_iq_wav(fn);
  x           = x(2:end); # omit the first chunk of 256 samples
  last_gpsfix = max(cat(1,x.gpslast));
  idx         = find(cat(1,x.gpslast) < max_last_gpsfix);
  xx          = {};
  if numel(idx) < 2
    return
  endif

  ## filter x
  x = x(idx);
  n = numel(x); ## from here on n is always >=2

  ## indices of fresh GNSS timestamps
  idx_ts = find(diff(cat(1, x.gpslast)) < 0);
  idx_ts = make_complete(idx_ts, n);

  ## remove outliers
  ts     = cat(1,x.gpssec) + 1e-9*cat(1,x.gpsnsec);
  idx_ts = fix_outliers(ts, idx_ts, fs);
  idx_ts = make_complete(idx_ts, n);
  ts     = ts(idx_ts);

  ## helper functions
  gpssec = @(s) s.gpssec + 1e-9*s.gpsnsec;
  get_fs = @(i) 512*(idx_ts(i) - idx_ts(i-1)) / (gpssec(x(idx_ts(i))) - gpssec(x(idx_ts(i-1))));

  ## keep history
  fs = [];

  ## (1) interval before 1st fresh timestamp: extrapolate backwards
  if idx_ts(1) > 1
    fs(end+1) = get_fs(2);
    k = 0;
    for j=1:idx_ts(1)-1
      xx(j).t = ts(1) - 512*(idx_ts(1)-1)/fs(end) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
      xx(j).z = x(j).z;
      k      += 1;
    end
  end

  ## (2) in between fresh timestamps: extrapolate between fresh timestamps
  for i=2:numel(idx_ts)
    ##__fs = [ idx_ts(i) - idx_ts(i-1) gpssec(x(idx_ts(i))) - gpssec(x(idx_ts(i-1))) ]
    fs(end+1) = get_fs(i);
    k = 0;
    for j=idx_ts(i-1):idx_ts(i)-1
      xx(j).t = ts(i-1) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
      xx(j).z = x(j).z;
      k      += 1;
    end
  end

  ## (3) after the last fresh timestamp: extrapolate after the last fresh timestamp
  k = 0;
  for j=idx_ts(end):n
    xx(j).t = ts(end) + (512*k+[0:numel(x(j).z)-1])'/fs(end);
    xx(j).z = x(j).z;
    k      += 1;
  end

  ## reduce the influence of outliers
  fs = median(fs);
endfunction

function idx_ts=fix_outliers(t, idx_ts, fs)
  ## more than one sample difference -> outlier
  outliers = find(abs(diff(t(idx_ts))*fs - 512*diff(idx_ts)) > 1);

  b = ones(size(idx_ts))==1;
  m = numel(outliers);
  for i=1:2:m
    if i+1 <= m
      b(outliers(i)+1:outliers(i+1)+1) = false;
    else
      b(outliers(i)+1:end) = false;
    end
  end
  ## filter idx_ts
  idx_ts(~b) = [];
endfunction

function idx_ts=make_complete(idx_ts, n)
  ## add fake fresh GNSS timestamps indices if numel(idx_ts)<2
  if isempty(idx_ts)
    idx_ts = [1 n];
  elseif numel(idx_ts) == 1
    if n-idx_ts > idx_ts
      idx_ts = [idx_ts n];
    else
      idx_ts = [1 idx_ts];
    end
  end
endfunction
