% ofdm_lib.m
% David Rowe Mar 2017

#{
  Library of functions that implement a PSK OFDM modem.
#}

1;
qam16;
esno_est;
ofdm_mode;
ofdm_state;
ofdm_helper;

%-------------------------------------------------------------
% ofdm_init
%-------------------------------------------------------------

#{

  Modem frame has a pilot every Ns symbols. There are Ns-1 data
  symbols between every pilot.

   e.g. for Ns=4, Nc=6:

    |-Nc-|          Time
    DDDDDD           |
   PPPPPPPP  ---     |
    DDDDDD    |      |
    DDDDDD    Ns     |
    DDDDDD    |      |
   PPPPPPPP  ---    \|/
    DDDDDD    |      |

   Freq------------------>

   In this figure, time flows down, freq across.
#}

function states = ofdm_init(config)
  Rs = config.Rs; 
  Tcp = config.Tcp; 
  Ns = config.Ns; 
  Nc = config.Nc;
  bps = config.bps;
  Np = config.Np;
  Ntxtbits = config.Ntxtbits;
  Nuwbits = config.Nuwbits;
  ftwindow_width = config.ftwindow_width;
  timing_mx_thresh = config.timing_mx_thresh;
  tx_uw = config.tx_uw;
  bad_uw_errors = config.bad_uw_errors;
  amp_scale = config.amp_scale;
  amp_est_mode = config.amp_est_mode;
  EsNo_est_all_symbols = config.EsNo_est_all_symbols;
  EsNodB = config.EsNodB; 
  state_machine = config.state_machine; 
  edge_pilots = config.edge_pilots;
  clip_gain1 = config.clip_gain1;
  clip_gain2 = config.clip_gain2;
  foff_limiter = config.foff_limiter; 
  txbpf_width_Hz = config.txbpf_width_Hz;
  data_mode = config.data_mode;
       
  states.Fs = 8000;
  states.bps = bps;
  states.Rs = Rs;
  states.Tcp = Tcp;
  states.Ns = Ns;                                 % one pilot every Ns symbols, e.g. Ns=4, ...PDDDPDDDP...
  states.Nc = Nc;                                 % Number of carriers
  states.M  = states.Fs/Rs;                       % oversampling rate
  states.Ncp = Tcp*states.Fs;
  states.Nbitsperframe = (Ns-1)*Nc*bps;           % total bits in all data symbols in modem frame
  states.Nsampersymbol = states.M+states.Ncp;     % number of samples in a single symbol
  states.Nsamperframe  = Ns*states.Nsampersymbol; % number of samples in a modem frame
  states.qam16 = [
    1 + j,  1 + j*3,  3 + j,  3 + j*3;
    1 - j,  1 - j*3,  3 - j,  3 - j*3;
   -1 + j, -1 + j*3, -3 + j, -3 + j*3;
   -1 - j, -1 - j*3, -3 - j, -3 - j*3]/3;
  rms = sqrt(states.qam16(:)'*states.qam16(:)/16);% set average Es to 1
  states.qam16 /= rms;
  states.qam16 *= exp(-j*pi/4);                   % same rotation as QPSK constellation
  states.Np = Np;                                 % number of modem frames per packet. In some modes we want
                                                  % the total packet of data to span multiple modem frames, e.g. HF data
                                                  % and/or when the FEC codeword is larger than the one
                                                  % modem frame.  In other modes (e.g. 700D/2020) Np=1, ie the modem frame
                                                  % is the same length as the packet/FEC codeword.
  states.Nbitsperpacket = Np*states.Nbitsperframe;
  states.Tpacket = Np*Ns*(Tcp+1/Rs);              % time for one packet in ms

  states.Ntxtbits = Ntxtbits;                     % reserved bits/frame for auxiliary text information.  Uncoded/unprotected so may
                                                  % be of limited use going forward, consider setting to 0
  states.Nuwbits  = Nuwbits;

  % some basic sanity checks
  assert(floor(states.M) == states.M);

  % UW symbol placement. 
  % Note we need to fill each UW symbols with bits.  The LDPC decoder
  % works on symbols so we can't break up any symbols into UW/FEC
  % encoded bits.

  states.uw_ind = states.uw_ind_sym = [];

  % lets see if all UW syms will fit in frame
  Nuwsyms = states.Nuwbits/bps;
  Ndatasymsperframe = (Ns-1)*Nc;
  states.spread_uw = 0;
  if states.spread_uw
    uw_step = 1.8*floor(states.Nbitsperpacket/states.Nuwbits);
  else
    uw_step = Nc+1;                                % default step for UW sym placement
  end  
  last_sym = floor(Nuwsyms*uw_step/bps+1);
  if last_sym > states.Np*Ndatasymsperframe
    uw_step = Nc-1;                                 % try a different step
  end
  last_sym = floor(Nuwsyms*uw_step/bps+1);
  assert(last_sym <= states.Np*Ndatasymsperframe);  % we still can't fit them all
  
  % Place UW symbols in frame
  for i=1:Nuwsyms
    ind_sym = floor(i*uw_step/bps+1);
    % printf("%d sym: %d\n",i, ind_sym);
    states.uw_ind_sym = [states.uw_ind_sym ind_sym];   % symbol index
    for b=bps-1:-1:0
      states.uw_ind = [states.uw_ind bps*ind_sym-b];   % bit index
    end
  end

  % how many of the first few frames have UW symbols in them
  Nsymsperframe = states.Nbitsperframe/states.bps;
  states.Nuwframes = ceil(states.uw_ind_sym(end)/Nsymsperframe);

  states.tx_uw = tx_uw;
  assert(length(states.tx_uw) == states.Nuwbits);
  tx_uw_syms = [];
  for b=1:bps:states.Nuwbits
    if bps == 2 tx_uw_syms = [tx_uw_syms qpsk_mod(states.tx_uw(b:b+1))]; end
    if bps == 4 tx_uw_syms = [tx_uw_syms qam16_mod(states.qam16, states.tx_uw(b:b+bps-1))]; end
  end
  states.tx_uw_syms = tx_uw_syms;
  
  % if the UW has this many errors it is "bad", the binomal cdf can be used to
  % set this with the ofdm_determine_bad_uw_errors() function below
  %
  %   Nuw=12; plot(0:Nuw, binocdf(0:Nuw,Nuw,0.05)); hold on; plot(binocdf(0:Nuw,Nuw,0.5)); hold off;
  states.bad_uw_errors = bad_uw_errors;

  states.ofdm_peak = 16384;
  % use this to scale tx output to 16 bit short to a peak value of 16384.  Adjusted by experiment
  states.amp_scale = amp_scale;
  % when using the clipping, this is the manual gain value.  Adjusted by experiment, trade off between
  % increased average power and BER
  states.clip_gain1 = clip_gain1;
  states.clip_gain2 = clip_gain2;
  states.txbpf_width_Hz = txbpf_width_Hz;

  % this is used to scale inputs to LDPC decoder to make it amplitude indep
  states.mean_amp = 0;

  % use a fixed EsNo for LDPC decoder, this seems to work OK and avoid another estimator
  states.EsNodB = EsNodB;

  % generate same BPSK pilots each time
  rand('seed',1);
  states.pilots = 1 - 2*(rand(1,Nc+2) > 0.5);
  %printf("number of pilots total: %d\n", length(states.pilots));

  % If set, place pilots at carrier 1 and Nc+2 to support low bandwidth phase est over grid
  % of 12 pilot_samples.  Used for 700D and 2020
  states.edge_pilots = edge_pilots;
  if states.edge_pilots == 0
     states.pilots(1) = 0;
     states.pilots(Nc+2) = 0;
  end
  
  % carrier tables for up and down conversion
  states.fcentre = fcentre = 1500;
  alower = fcentre - Rs * (Nc/2);  % approx frequency of lowest carrier
  Nlower = round(alower / Rs) - 1; % round this to nearest integer multiple from 0Hz to keep DFT happy
  %printf("  fcentre: %f alower: %f alower/Rs: %f Nlower: %d\n", fcentre, alower, alower/Rs, Nlower);
  w = (Nlower:Nlower+Nc+1)*2*pi/(states.Fs/Rs);
  W = zeros(Nc+2,states.M);
  for c=1:Nc+2
    W(c,:) = exp(j*w(c)*(0:states.M-1));
  end
  states.w = w;
  states.W = W;

  % fine timing search +/- window_width/2 from current timing instant,
  % set this to roughly twice the maximum delay spread
  states.ftwindow_width = ftwindow_width;

  % magic number we adjust by experiment (see ofdm_dev.m acquisition tests, blog post on 700D sync)
  states.timing_mx_thresh = timing_mx_thresh;

  % Receive buffer: rxbufst + D P DDD P DDD P DDD P D
  %                                   ^
  %                                   nominal start of current modem frame
  
  if length(data_mode)
    Nrxbufhistory = (states.Np+2)*states.Nsamperframe;     % extra storage at start of rxbuf to allow us to step back in time 
  else
    Nrxbufhistory = 0;
  end
  states.rxbufst = Nrxbufhistory;                          % start of rxbuf window used for demod of current rx frame
  states.Nrxbufhistory = Nrxbufhistory;
  
  %                       D                    P DDD P DDD P DDD             P                    D
  states.Nrxbufmin = states.Nsampersymbol + 3*states.Nsamperframe + states.Nsampersymbol + states.Nsampersymbol;
  states.Nrxbuf = Nrxbufhistory + states.Nrxbufmin;
  states.rxbuf = zeros(1, states.Nrxbuf);

  % default settings on a bunch of options and states

  states.verbose = 0;
  states.timing_en = 1;
  states.foff_est_en = 1;
  states.phase_est_en = 1;
  states.phase_est_bandwidth = "high";
  states.dpsk = 0;
  states.amp_est_mode = amp_est_mode;

  states.foff_est_gain = 0.1;
  states.foff_limiter = foff_limiter;
  states.foff_est_hz = 0;
  states.sample_point = states.timing_est = 1;
  states.nin = states.Nsamperframe;
  states.timing_valid = 0;
  states.timing_mx = 0;
  states.coarse_foff_est_hz = 0;

  states.foff_metric = 0;

  % generate OFDM pilot symbol, used for timing and freq offset est

  rate_fs_pilot_samples = states.pilots * W/states.M;

  % During tuning it was found that not including the cyc prefix in
  % rate_fs_pilot_samples produced better fest results

  %states.rate_fs_pilot_samples = [rate_fs_pilot_samples(states.M-states.Ncp+1:states.M) rate_fs_pilot_samples];
  states.rate_fs_pilot_samples = [zeros(1,states.Ncp) rate_fs_pilot_samples];

  % pre-compute a constant used to detect valid modem frames

  Npsam = length(states.rate_fs_pilot_samples);
  states.timing_norm = Npsam*(states.rate_fs_pilot_samples * states.rate_fs_pilot_samples');
  % printf("timing_norm: %f\n", states.timing_norm)

  % sync state machine

  states.sync_state = states.last_sync_state = 'search';
  states.uw_errors = 0;
  states.sync_counter = 0;
  states.frame_count = 0;                                 % number of frames we have been in sync
  states.sync_start = 0;
  states.sync_end = 0;
  states.modem_frame = 0;                                 % keep track of how many frames received in packet
  states.state_machine = state_machine;                   % mode specific state machine
  states.packetsperburst = 0;                             % for OFDM data modes, how many packets before we reset state machine
  states.postambledetectoren = strcmp(data_mode,"burst");
  states.npre = states.npost = 0;                         % counters for logging
  
  % LDPC code is optionally enabled

  states.rate = 1.0;
  states.ldpc_en = 0;

  % init some output states for logging

  states.rx_sym = zeros(1+Ns+1+1, Nc+2);

  % Es/No (SNR) est states

  states.EsNo_est_all_symbols = EsNo_est_all_symbols;
  states.clock_offset_est = 0;

  % pre-amble for data modes
  states.data_mode = data_mode;
  if length(states.data_mode)
    states.tx_preamble = ofdm_generate_preamble(states, 2);
    states.tx_postamble = ofdm_generate_preamble(states, 3);
  end
  
  % automated tests
  test_qam16_mod_demod(states.qam16);
  test_assemble_disassemble(states);
endfunction

% Gray coded QPSK modulation function
function symbol = qpsk_mod(two_bits)
    two_bits_decimal = sum(two_bits .* [2 1]);
    switch(two_bits_decimal)
        case (0) symbol =  1;
        case (1) symbol =  j;
        case (2) symbol = -j;
        case (3) symbol = -1;
    endswitch
endfunction

% Gray coded QPSK demodulation function
function two_bits = qpsk_demod(symbol)
    bit0 = real(symbol*exp(j*pi/4)) < 0;
    bit1 = imag(symbol*exp(j*pi/4)) < 0;
    two_bits = [bit1 bit0];
endfunction


function out = freq_shift(in, foff, Fs)
  foff_rect = exp(j*2*pi*foff/Fs);
  foff_phase_rect = exp(j*0);

  for r=1:length(in)
    foff_phase_rect *= foff_rect;
    out(r) = in(r)*foff_phase_rect;
  end
endfunction


% -----------------------------------------------------------------
% ofdm_mod - modulates a complete packet (one or more modem frames)
% ----------------------------------------------------------------

function tx = ofdm_mod(states, tx_bits)
  ofdm_load_const;
  assert(length(tx_bits) == Nbitsperpacket);

  % map to symbols in linear array

  if bps == 1
    tx_sym_lin = 2*tx_bits - 1;
  end
  if bps == 2
    for s=1:Nbitsperpacket/bps
      tx_sym_lin(s) = qpsk_mod(tx_bits(2*(s-1)+1:2*s));
    end
  end
  if bps == 4
    for s=1:Nbitsperpacket/bps
      tx_sym_lin(s) = qam16_mod(states.qam16,tx_bits(4*(s-1)+1:4*s));
    end
  end

  tx = ofdm_txframe(states, tx_sym_lin);
endfunction


% ----------------------------------------------
% ofdm_txframe - modulates one packet of symbols
% ----------------------------------------------

function tx = ofdm_txframe(states, tx_sym_lin)
  ofdm_load_const;
  assert(length(tx_sym_lin) == Nbitsperpacket/bps);

  % place data symbols in multi-carrier frame with pilots and boundary carriers

  s = 1; tx_frame = zeros(Np*Ns,Nc+2);
  for r=1:Np*Ns
    if mod(r-1,Ns) == 0
      % row of pilots
      tx_frame(r,:) = pilots;
    else
      % row of data symbols
      arowofsymbols = tx_sym_lin(s:s+Nc-1);
      tx_frame(r,2:Nc+1) = arowofsymbols;
      s += Nc;
      if states.dpsk
        tx_frame(r,2:Nc+1) = tx_frame(r,2:Nc+1) .* tx_frame(r-1,2:Nc+1);
      end
    end
  end
  % make sure we use all the symbols
  assert((s-1) == length(tx_sym_lin));

  % OFDM upconvert symbol by symbol so we can add CP

  tx = [];
  for r=1:Ns*Np
    asymbol = tx_frame(r,:) * W/M;
    asymbol_cp = [asymbol(M-Ncp+1:M) asymbol];
    tx = [tx asymbol_cp];
  end
endfunction


% -----------------------------------------------------------
% est_timing
% -----------------------------------------------------------

#{
  Correlates known samples (for example pilots or a preamble) with a window of received
  samples to determine the most likely timing offset.  Optionally combines
  known samples from two frames (e.g. pilots at start of this and next frame)
  so we need at least Nsamperframe+M+Ncp samples in rx.

  Can be used for acquisition (coarse timing), and fine timing.  Tends
  to break down when freq offset approaches +/- symbol rate (e.g +/-
  25 Hz for 700D).
#}

function [t_est timing_valid timing_mx av_level] = est_timing(states, rx, known_samples, step, dual=1)
    ofdm_load_const;
    Npsam = length(known_samples);

    Ncorr = length(rx) - (Nsamperframe+Npsam);
    corr = zeros(1,Ncorr);
    %printf("Npsam: %d M+Ncp: %d Ncorr: %d Nsamperframe: %d step: %d\n", Npsam,  M+Ncp, Ncorr, Nsamperframe, step);

    % normalise correlation so we can compare to a threshold across varying input levels

    av_level = 2*sqrt(states.timing_norm*(rx*rx')/length(rx)) + 1E-12;

    % correlate with pilots at start and (optionally) end of frame to determine timing offset

    for i=1:step:Ncorr
      rx1 = rx(i:i+Npsam-1); 
      corr_st = rx1 * known_samples'; 
      corr_en = 0;
      if dual
        % for the streaming voice modes we also correlate with pilot samples at start of next frame 
        rx2 = rx(i+Nsamperframe:i+Nsamperframe+Npsam-1);
        corr_en = rx2 * known_samples';
      end
      corr(i) = (abs(corr_st) + abs(corr_en))/av_level;
    end

    [timing_mx t_est] = max(abs(corr));
    % only declare timing valid if there are enough samples in rxbuf to demodulate a frame
    timing_valid = (abs(rx(t_est)) > 0) && (timing_mx > timing_mx_thresh);
    
    if verbose > 1
      printf("  av_level: %5.4f mx: %4.3f timing_est: %4d timing_valid: %d\n", av_level, timing_mx, t_est, timing_valid);
    end
    if verbose > 2
      figure(10); clf;
      subplot(211); plot(rx)
      subplot(212); plot(corr)
      figure(11); clf; plot(real(known_samples));
    end

endfunction


% -----------------------------------------------------------
% est_freq_offset_known_corr
% -----------------------------------------------------------

#{
  Determines frequency offset at current timing estimate, used for
  coarse freq offset estimation during streaming mode acquisition.
#}

function foff_est = est_freq_offset_known_corr(states, rx, known_samples, t_est, dual=1)
    ofdm_load_const;
    Npsam = length(known_samples);

    % extract pilot samples from either end of frame
    rx1  = rx(t_est:t_est+Npsam-1); rx2 = rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam-1);

    % "mix" these down (correlate) with 0 Hz offset pilot samples
    corr_st = rx1 .* conj(known_samples);
    if dual 
      corr_en = rx2 .* conj(known_samples);
    end
    
    % sample sum of DFT magnitude of correlated signals at each freq offset and look for peak
    st = -20; en = 20; foff_est = 0; Cabs_max = 0;

    for f=st:en
       w = 2*pi*f/Fs;
       C_st = corr_st * exp(j*w*(0:Npsam-1))';
       C_en = 0;
       if dual
          C_en = corr_en * exp(j*w*(0:Npsam-1))';
       end
       Cabs = abs(C_st) + abs(C_en);
       %printf("f: %4.1f Cabs: %f Cmax: %f\n", f, Cabs, Cabs_max);
       if Cabs > Cabs_max
         Cabs_max = Cabs;
         foff_est = f;
       end
    end

    if states.verbose > 1
      printf("  foff_est: %f\n", foff_est);
    end

endfunction


% Joint estimation used for data mode burst acquistion

function [t_est foff_est timing_mx] = est_timing_and_freq(states, rx, known_samples, tstep, fmin, fmax, fstep)
    ofdm_load_const;
    Npsam = length(known_samples);

    Ncorr = length(rx) - Npsam + 1;
    corr = zeros(1,Ncorr);
    
    % set up matrix of freq shifted known samples for correlation with received signal.  Each row
    % is the known samples shifted by a different freq offset
    
    M = [];
    for afcoarse=fmin:fstep:fmax
       w = 2*pi*afcoarse/Fs;
       wvec = exp(j*w*(0:Npsam-1));
       M = [M; known_samples .* wvec];
    end
    
    % At each timing position, correlate with known samples at all possible freq offsets.  Result
    % is a column vector for each timing offset.  Each matrix cell is a freq,timing coordinate
    
    corr = [];
    for t=1:tstep:Ncorr
      rx1 = rx(t:t+Npsam-1); 
      col = M * rx1';
      corr = [corr, col];
    end
    
    % best timing offset is the col with the global max of the corr matrix
    max_col = max(abs(corr));
    [mx mx_col] = max(max_col);
    t_est = (mx_col-1)*tstep;
    
    % obtain normalised real number for timing mx
    mag1 = known_samples*known_samples';
    mag2 = rx(t_est+1:t_est+Npsam)*rx(t_est+1:t_est+Npsam)';
    timing_mx = mx*mx'/(mag1*mag2+1E-12);
    
    % determine frequency offset for row where max is located
    [tmp freq_row] = max(corr(:,mx_col));
    foff_est = fmin + fstep*(freq_row-1);
       
    if verbose > 1
      printf("  t_est: %d timing:mx: %f foff_est: %f\n", t_est, timing_mx, foff_est);
    end
    if verbose > 2
      figure(10); clf;
      subplot(211); plot(rx)
      subplot(212); plot(corr)
      figure(11); clf; plot(real(known_samples));
    end

endfunction


% streaming mode acquistion, used mainly for voice modes

function [timing_valid states] = ofdm_sync_search_stream(states)
    ofdm_load_const;
    
    st = rxbufst + M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe + M+Ncp - 1;
  
    % Attempt coarse timing estimate (i.e. detect start of frame) at a range of frequency offsets

    timing_mx = 0; fcoarse = 0; timing_valid = 0; ct_est = 1;
    for afcoarse=-40:40:40
      % vector of local oscillator samples to shift input vector
      % these could be computed on the fly to save memory, or pre-computed in flash at tables as they are static

      if afcoarse != 0
        w = 2*pi*afcoarse/Fs;
        wvec = exp(-j*w*(0:2*Nsamperframe+M+Ncp-1));

        % choose best timing offset metric at this freq offset
        [act_est atiming_valid atiming_mx] = est_timing(states, wvec .* states.rxbuf(st:en), states.rate_fs_pilot_samples, 2);
      else
        % exp(-j*0) is just 1 when afcoarse is 0
        [act_est atiming_valid atiming_mx] = est_timing(states, states.rxbuf(st:en), states.rate_fs_pilot_samples, 2);
      end

      %printf("afcoarse: %f atiming_mx: %f\n", afcoarse, atiming_mx);

      if atiming_mx > timing_mx
        ct_est = act_est;
        timing_valid = atiming_valid;
        timing_mx = atiming_mx;
        fcoarse = afcoarse;
      end
    end

    % refine freq est within -/+ 20 Hz window

    if fcoarse != 0
      w = 2*pi*fcoarse/Fs;
      wvec = exp(-j*w*(0:2*Nsamperframe+M+Ncp-1));
      foff_est = est_freq_offset_known_corr(states, wvec .* states.rxbuf(st:en), states.rate_fs_pilot_samples, ct_est);
      foff_est += fcoarse;
    else
      % exp(-j*0) is just 1 when fcoarse is 0
      foff_est = est_freq_offset_known_corr(states, states.rxbuf(st:en), states.rate_fs_pilot_samples, ct_est);
    end

    if verbose
      printf(" ct_est: %4d mx: %3.2f coarse_foff: %5.1f timing_valid: %d", ct_est, timing_mx, foff_est, timing_valid);
    end
    
  if timing_valid
    states.nin = ct_est - 1;
  else
    states.nin = Nsamperframe;
  end

  states.timing_valid = timing_valid;
  states.timing_mx = timing_mx;
  states.coarse_foff_est_hz = foff_est;
  states.sample_point = states.timing_est = 1;
endfunction
 

% two stage acquisition detector for burst mode

function results = burst_acquisition_detector(states, rx, n, known_sequence)
  ofdm_load_const;
    
  % initial search over coarse grid
  tstep = 4; fstep = 5;
  [ct_est foff_est timing_mx] = est_timing_and_freq(states, rx(n:n+2*Nsamperframe-1), known_sequence, 
                                                    tstep, fmin = -50, fmax = 50, fstep);
  % refine estimate over finer grid                             
  fmin = foff_est - ceil(fstep/2); fmax = foff_est + ceil(fstep/2); 
  fine_st = max(1, n + ct_est - tstep/2); fine_en = fine_st + Nsamperframe + tstep - 1;
  [ct_est foff_est timing_mx] = est_timing_and_freq(states, rx(fine_st:fine_en), known_sequence, 1, fmin, fmax, 1);
  % refer ct_est to nominal start of frame rx_buf(n)
  ct_est += fine_st - n;
  results.ct_est = ct_est; results.foff_est = foff_est; results.timing_mx = timing_mx;
end


% Burst mode acquisition ------------------------------------------

function [timing_valid states] = ofdm_sync_search_burst(states)
  ofdm_load_const;

  pre_post = "";
  st = rxbufst + M+Ncp + Nsamperframe + 1; en = st + 2*Nsamperframe - 1;
  pre = burst_acquisition_detector(states, states.rxbuf, st, states.tx_preamble);
  if states.postambledetectoren
    post = burst_acquisition_detector(states, states.rxbuf, st, states.tx_postamble);
  end
  
  if isfield(states,"postambletest") pre.timing_mx = 0; end % force ignore preamble to test postamble

  if (states.postambledetectoren == 0) || (pre.timing_mx > post.timing_mx)
    timing_mx = pre.timing_mx; ct_est = pre.ct_est; foff_est = pre.foff_est;
    pre_post = "pre";
  else
    timing_mx = post.timing_mx; ct_est = post.ct_est; foff_est = post.foff_est;
    pre_post = "post";
  end
  timing_valid = timing_mx > timing_mx_thresh;
  
  if timing_valid
    % potential candidate found ....

    % calculate number of samples we need on next buffer to get into sync
    if strcmp(pre_post, "post")
      states.nin = 0;
      % printf("\n  rxbufst: %d ", states.rxbufst);
      states.rxbufst -= states.Np*states.Nsamperframe; % backup to first modem frame in packet
      states.rxbufst += ct_est - 1;
      states.npost++;
      % printf("%d\n", states.rxbufst);
    else
      % ct_est is start of preamble, so advance past that to start of first modem frame
      states.nin = Nsamperframe + ct_est - 1;
      states.npre++;
    end
  else
    states.nin = Nsamperframe;
  end

  states.ct_est = ct_est;
  states.timing_valid = timing_valid;
  states.timing_mx = timing_mx;
  states.sample_point = states.timing_est = 1;
  states.foff_est_hz = foff_est;

  if verbose
    printf("  ct_est: %4d nin: %4d mx: %3.2f foff_est: %5.1f timing_valid: %d %4s", 
           ct_est, states.nin, timing_mx, foff_est, timing_valid, pre_post);
  end
endfunction


% ----------------------------------------------------------------------------------
% ofdm_sync_search - attempts to find coarse sync parameters for modem initial sync
% ----------------------------------------------------------------------------------

function [timing_valid states] = ofdm_sync_search(states, rxbuf_in)
  ofdm_load_const;

  % update rxbuf so it is primed for when we have to call ofdm_demod()

  states.rxbuf(1:Nrxbuf-states.nin) = states.rxbuf(states.nin+1:Nrxbuf);
  states.rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in;
  
  if strcmp(states.data_mode, "burst")
    [timing_valid states] = ofdm_sync_search_burst(states);
  else
    [timing_valid states] = ofdm_sync_search_stream(states);
  end
endfunction


% ------------------------------------------
% ofdm_demod - Demodulates one frame of bits
% ------------------------------------------

#{

  For phase estimation we need to maintain buffer of 3 frames plus
  one pilot, so we have 4 pilots total. '^' is the start of current
  frame that we are demodulating.

  P DDD P DDD P DDD P
        ^

  Then add one symbol either side to account for movement in
  sampling instant due to sample clock differences:

  D P DDD P DDD P DDD P D
          ^

  Returns:
    rx_bits    - (hard decoded/raw/uncoded) demodulated data bits from packet
    aphase_est - phase est for each data symbol
    rx_np      - output data symbols after phase correction
    rx_amp     - amplitude estimates for each symbol
#}

function [states rx_bits achannel_est_rect_log rx_np rx_amp] = ofdm_demod(states, rxbuf_in)
  ofdm_load_const;

  % insert latest input samples into rxbuf

  rxbuf(1:Nrxbuf-states.nin) = rxbuf(states.nin+1:Nrxbuf);
  rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in;

  % get latest freq offset estimate

  woff_est = 2*pi*foff_est_hz/Fs;

  % update timing estimate --------------------------------------------------

  delta_t = coarse_foff_est_hz = timing_valid = timing_mx = 0;
  if timing_en
    % update timing at start of every frame

    % search for timing in a window centered on timing_est, the window will typically be around 2Ncp wide as we could
    % get a shift of +Ncp or -Ncp if we swing from one delay extreme to another
    st = rxbufst + M+Ncp + Nsamperframe + 1 - floor(ftwindow_width/2) + (timing_est-1);
    en = st + Nsamperframe-1 + M+Ncp + ftwindow_width-1;

    [ft_est timing_valid timing_mx] = est_timing(states, rxbuf(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples, 1);
    % printf("  timing_est: %d ft_est: %d timing_valid: %d timing_mx: %d\n", timing_est, ft_est, timing_valid, timing_mx);

    % if we are in a deep fade timing_valid will not be asserted as ft_est will be garbage, so we don't
    % adjust timing est, just freewheel for now
    if timing_valid
        
      % adjust timing_est based on ft_est    
      timing_est = timing_est + ft_est - ceil(ftwindow_width/2);

      % Track the ideal sampling point, which is Ncp for a multipath signal whose delay varies between 0 and Ncp.  The
      % timing est will be bouncing back and forth due to multipath so we may need to use the upper or lower limit of
      % the timing est to track the ideal sample_point. A good way to explore this algorithm is to disable the feedback
      % loop for nin adjustment below, and look at the plots from ofdm_rx with +ve and -ve sample clock offsets 
      % (sox can be used to resample).  The "4" constants are small guard bands so we don't stumble outside of the CP 
      % due to noise.
      
      delta_t = ft_est - ceil(ftwindow_width/2);           % just used for plotting
      sample_point = max(timing_est+4, sample_point);      % we are at max timing est, so sample point just above
      sample_point = min(timing_est+Ncp-4, sample_point);  % we are at min timing_est, so sample point Ncp above
    end

    if verbose > 1
      printf("  ft_est: %2d mx: %3.2f coarse_foff: %4.1f foff: %4.1f\n", ft_est, timing_mx, coarse_foff_est_hz, foff_est_hz);
    end

  end

  % down convert at current timing instant----------------------------------

   rx_sym = zeros(1+Ns+1+1, Nc+2);

  % previous pilot

  st = rxbufst + M+Ncp + Nsamperframe + (-Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1;

  for c=1:Nc+2
    acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
    rx_sym(1,c) = sum(acarrier);
  end

  % pilot - this frame - pilot

  for rr=1:Ns+1
    st = rxbufst + M+Ncp + Nsamperframe + (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1;
    for c=1:Nc+2
      acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
      rx_sym(rr+1,c) = sum(acarrier);
    end
  end

  % next pilot

  st = rxbufst + M+Ncp + Nsamperframe + (2*Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1;
  for c=1:Nc+2
    acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
    rx_sym(Ns+3,c) = sum(acarrier);
  end

  % est freq err based on all carriers ------------------------------------

  if foff_est_en
    freq_err_rect = sum(rx_sym(2,:))' * sum(rx_sym(2+Ns,:));

    % prevent instability in atan(im/re) when real part near 0

    freq_err_rect += 1E-6;

    %printf("freq_err_rect: %f %f angle: %f\n", real(freq_err_rect), imag(freq_err_rect), angle(freq_err_rect));
    freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns);
    if states.foff_limiter
      freq_err_hz = max(freq_err_hz,-1);
      freq_err_hz = min(freq_err_hz, 1);
    end
    foff_est_hz = foff_est_hz + foff_est_gain*freq_err_hz;
  end

  % OK - now channel for each carrier and correct phase  ----------------------------------

  achannel_est_rect = zeros(1,Nc+2);
  aamp_est_pilot = zeros(1,Nc+2);
  for c=2:Nc+1

    % estimate channel for this carrier using an average of 12 pilots
    % in a rect 2D window centred on this carrier

    % PPP  <-- frame-1
    % ---
    % PPP  <-- you are here
    % DDD
    % DDD
    % PPP  <-- frame+1
    % ---
    % PPP  <-- frame+2

    if isfield(states, "phase_est_bandwidth")
      phase_est_bandwidth = states.phase_est_bandwidth;
    else
      phase_est_bandwidth = "low";
    end

    if strcmp(phase_est_bandwidth, "high")
      % Only use pilots at start and end of this frame to track quickly changes in phase
      % present.  Useful for initial sync where freq offset est may be a bit off, and
      % for high Doppler channels.  As less pilots are averaged, low SNR performance
      % will be poorer.
      achannel_est_rect(c) =  rx_sym(2,c)*pilots(c)';        % frame
      achannel_est_rect(c) += rx_sym(2+Ns,c)*pilots(c)';     % frame+1
      aamp_est_pilot(c) = abs(rx_sym(2,c)) + abs(rx_sym(2+Ns,c));
    elseif strcmp(phase_est_bandwidth, "low")
      % Average over a bunch of pilots in adjacent carriers, and past and future frames, good
      % low SNR performance, but will fall over with high Doppler or freq offset.
      cr = c-1:c+1;
      achannel_est_rect(c) =  rx_sym(2,cr)*pilots(cr)';      % frame
      achannel_est_rect(c) += rx_sym(2+Ns,cr)*pilots(cr)';   % frame+1
      aamp_est_pilot(c)  = sum(abs(rx_sym(2,cr)));
      aamp_est_pilot(c) += sum(abs(rx_sym(2+Ns,cr)));

      % use next step of pilots in past and future

      achannel_est_rect(c) += rx_sym(1,cr)*pilots(cr)';      % frame-1
      achannel_est_rect(c) += rx_sym(2+Ns+1,cr)*pilots(cr)'; % frame+2
      aamp_est_pilot(c) += sum(abs(rx_sym(1,cr)));
      aamp_est_pilot(c) += sum(abs(rx_sym(2+Ns+1,cr)));
    end
  end

  % pilots are estimated over 12 pilot symbols, so find average

  if strcmp(phase_est_bandwidth, "high")
    achannel_est_rect /= 2;
    aamp_est_pilot /= 2;
  elseif strcmp(phase_est_bandwidth, "low")
    achannel_est_rect /= 12;
    aamp_est_pilot /= 12;
  end

  aphase_est_pilot = angle(achannel_est_rect);
  if states.amp_est_mode == 0
    % legacy 700D/2020 ampl estimator for compatibility with current C code
    aamp_est_pilot = abs(achannel_est_rect);
  end
  achannel_est_rect = aamp_est_pilot.*exp(j*aphase_est_pilot);

  % correct phase offset using phase estimate, and demodulate
  % bits, separate loop as it runs across cols (carriers) to get
  % frame bit ordering correct

  rx_bits = []; rx_np = []; rx_amp = []; achannel_est_rect_log = [];
  for rr=1:Ns-1
    for c=2:Nc+1
      if phase_est_en
        if states.dpsk
          rx_corr = rx_sym(rr+2,c) *  rx_sym(rr+1,c)';
        else
          rx_corr = rx_sym(rr+2,c) * exp(-j*aphase_est_pilot(c));
        end
      else
        rx_corr = rx_sym(rr+2,c);
      end

      rx_np = [rx_np rx_corr];
      rx_amp = [rx_amp aamp_est_pilot(c)];

      % hard decision demod
      if bps == 1 abit = real(rx_corr) > 0; end
      if bps == 2 abit = qpsk_demod(rx_corr); end
      if bps == 4 abit = qam16_demod(states.qam16, rx_corr, max(1E-12,aamp_est_pilot(c))); end
      rx_bits = [rx_bits abit];
    end % c=2:Nc+1
    achannel_est_rect_log = [achannel_est_rect_log; achannel_est_rect(2:Nc+1)];
  end

  % Adjust nin to take care of sample clock offset.  When debugong or exploring how timing loop works
  % it's a good idea to comment out this code to "open the loop".

  nin = Nsamperframe;
  if timing_en && timing_valid
    states.clock_offset_est = 0.9*states.clock_offset_est + 0.1*abs(states.timing_est - timing_est)/Nsamperframe;
    thresh = (M+Ncp)/8;
    tshift = (M+Ncp)/4;
    if timing_est > thresh
      nin = Nsamperframe+tshift;
      timing_est -= tshift;
      sample_point -= tshift;
    end
    if timing_est < -thresh
      nin = Nsamperframe-tshift;
      timing_est += tshift;
      sample_point += tshift;
    end
  end

  % use internal rxbuf samples if they are available
  rxbufst_next = rxbufst + nin;
  %printf("\nrxbufst: %d rxbufst_next: %d nin: %d Nrxbufmin: %d rqd: %d Nrxbuf: %d\n", 
  %     rxbufst, rxbufst_next, nin, Nrxbufmin, rxbufst_next + Nrxbufmin, Nrxbuf);
  if rxbufst_next + Nrxbufmin <= Nrxbuf
     % printf("Can maybe use rxbufst!\n");
     rxbufst = rxbufst_next;
     nin = 0;
  end
      
  % maintain mean amp estimate for LDPC decoder
  states.mean_amp = 0.9*states.mean_amp + 0.1*mean(rx_amp);

  states.rx_sym = rx_sym;
  states.rxbuf = rxbuf;
  states.nin = nin;
  states.rxbufst = rxbufst;
  states.timing_valid = timing_valid;
  states.timing_mx = timing_mx;
  states.timing_est = timing_est;
  states.sample_point = sample_point;
  states.delta_t = delta_t;
  states.foff_est_hz = foff_est_hz;
  states.coarse_foff_est_hz = coarse_foff_est_hz; % just used for tofdm
endfunction


function SNR3kdB = snr_from_esno(states, EsNodB)
    ofdm_load_const;

    % We integrate over M samples to get the received symbols.  Additional signal power
    % is used for the cyclic prefix samples.
    cyclic_power = 10*log10((Ncp+M)/M);
    % Es is the energy for each symbol.  To get signal power lets
    % multiply by symbols/second, and calculate noise power in 3000 Hz.
    SNR3kdB = EsNodB + 10*log10(Nc*Rs/3000) + cyclic_power;
endfunction

% ----------------------------------------------------------------------------------
% assemble_modem_packet - assemble modem packet from UW, payload, and txt bits
% ----------------------------------------------------------------------------------

function modem_frame = assemble_modem_packet(states, payload_bits, txt_bits)
  ofdm_load_const;

  # Due to the operation of the FEC encoder or interleaver, Tx data
  # usually comes in "packet size" chunks, so assembly operates on an
  # entire packet (multiple modem frames if Np>1)

  p = 1; u = 1;
  modem_frame = zeros(1,Nbitsperpacket);

  for b=1:Nbitsperpacket-Ntxtbits;
    if (u <= Nuwbits) && (b == uw_ind(u))
      modem_frame(b) = tx_uw(u++);
    else
      modem_frame(b) = payload_bits(p++);
    end
  end
  t = 1;
  for b=Nbitsperpacket-Ntxtbits+1:Nbitsperpacket
    modem_frame(b) = txt_bits(t++);
  end
  assert(u == (Nuwbits+1));
  assert(p = (length(payload_bits)+1));
endfunction


% ----------------------------------------------------------------------------------
% assemble_modem_packet_symbols - assemble modem packet from UW, payload, and txt bits
% ----------------------------------------------------------------------------------

function modem_frame = assemble_modem_packet_symbols(states, payload_syms, txt_syms)
  ofdm_load_const;

  Nsymsperpacket = Nbitsperpacket/bps;
  Nuwsyms = Nuwbits/bps;
  Ntxtsyms = Ntxtbits/bps;
  modem_frame = zeros(1,Nsymsperpacket);
  p = 1; u = 1;

  for s=1:Nsymsperpacket-Ntxtsyms;
    if (u <= Nuwsyms) && (s == uw_ind_sym(u))
      modem_frame(s) = states.tx_uw_syms(u++);
    else
      modem_frame(s) = payload_syms(p++);
    end
  end
  t = 1;
  for s=Nsymsperpacket-Ntxtsyms+1:Nsymsperpacket
    modem_frame(s) = txt_syms(t++);
  end
  assert(u == (Nuwsyms+1));
  assert(p = (length(payload_syms)+1));
endfunction


% ------------------------------------------------------------------------------------------------
% extract_uw - extract just the UW from the first few frames of a packet, to check UW
%              during acquisition
% -------------------------------------------------------------------------------------------------

function rx_uw = extract_uw(states, rx_syms, rx_amps)
  ofdm_load_const;

  Nsymsperframe = Nbitsperframe/bps;
  assert(length(rx_syms) == Nuwframes*Nsymsperframe);
  Nuwsyms = Nuwbits/bps;
  rx_uw_syms = zeros(1,Nuwsyms);
  rx_uw_amps = zeros(1,Nuwsyms);
  u = 1;

  for s=1:Nuwframes*Nsymsperframe
    if (u <= Nuwsyms) && (s == uw_ind_sym(u))
      rx_uw_syms(u) = rx_syms(s);
      rx_uw_amps(u) = rx_amps(s);
      u++;
    end
  end
  assert(u == (Nuwsyms+1));

  % now demodulate UW bits
  rx_uw = zeros(1,Nuwbits);

  for s=1:Nuwsyms
    if bps == 2
      rx_uw(bps*(s-1)+1:bps*s) = qpsk_demod(rx_uw_syms(s));
    elseif bps == 4
      rx_uw(bps*(s-1)+1:bps*s) = qam16_demod(states.qam16,rx_uw_syms(s), max(1E-12,rx_amps(s)));
    end
  end
endfunction


% ------------------------------------------------------------------------------------------------
% disassemble_modem_packet - extract UW, txt bits, and payload symbols from a packet of symbols
% -------------------------------------------------------------------------------------------------

function [rx_uw payload_syms payload_amps txt_bits] = disassemble_modem_packet(states, modem_frame_syms, modem_frame_amps)
  ofdm_load_const;

  Nsymsperpacket = Nbitsperpacket/bps;
  Nuwsyms = Nuwbits/bps;
  Ntxtsyms = Ntxtbits/bps;
  payload_syms = zeros(1,Nsymsperpacket-Nuwsyms-Ntxtsyms);
  payload_amps = zeros(1,Nsymsperpacket-Nuwsyms-Ntxtsyms);
  rx_uw_syms = zeros(1,Nuwsyms);
  rx_uw_amps = zeros(1,Nuwsyms);
  txt_syms = zeros(1,Ntxtsyms);
  p = 1; u = 1;

  for s=1:Nsymsperpacket-Ntxtsyms;
    if (u <= Nuwsyms) && (s == uw_ind_sym(u))
      rx_uw_syms(u) = modem_frame_syms(s);
      rx_uw_amps(u) = modem_frame_amps(s);
      u++;
    else
      payload_syms(p) = modem_frame_syms(s);
      payload_amps(p++) = modem_frame_amps(s);
    end
  end
  t = 1;
  for s=Nsymsperpacket-Ntxtsyms+1:Nsymsperpacket
    txt_syms(t++) = modem_frame_syms(s);
  end
  assert(u == (Nuwsyms+1));
  assert(p = (Nsymsperpacket+1));

  % now demodulate UW and txt bits

  rx_uw = zeros(1,Nuwbits);
  txt_bits = zeros(1,Ntxtbits);

  for s=1:Nuwsyms
    if bps == 2
      rx_uw(bps*(s-1)+1:bps*s) = qpsk_demod(rx_uw_syms(s));
    elseif bps == 4
      rx_uw(bps*(s-1)+1:bps*s) = qam16_demod(states.qam16,rx_uw_syms(s),rx_uw_amps(s));
    end
  end
  for s=1:Ntxtsyms
    txt_bits(2*s-1:2*s) = qpsk_demod(txt_syms(s));
  end

endfunction


%-----------------------------------------------------------------------
% ofdm_rand - a psuedo-random number generator that we can implement
%             in C with identical results to Octave.  Returns an unsigned
%             int between 0 and 32767
%-----------------------------------------------------------------------

function r = ofdm_rand(n, seed=1)
  r = zeros(1,n);
  for i=1:n
    seed = mod(1103515245 * seed + 12345, 32768);
    r(i) = seed;
  end
endfunction


% build a single modem frame preamble vector for reliable single frame acquisition
% on data modes
function tx_preamble = ofdm_generate_preamble(states, seed=2)
  tmp_states = states;
  % tweak local copy of states so we can generate a 1 modem-frame packet
  tmp_states.Np = 1; tmp_states.Nbitsperpacket = tmp_states.Nbitsperframe;
  preamble_bits = ofdm_rand(tmp_states.Nbitsperframe, seed) > 16384;
  tx_preamble = ofdm_mod(tmp_states, preamble_bits);
endfunction


% ------------------------------------------------------------------------------
% Handle FEC encoding/decoding
% ------------------------------------------------------------------------------

function [frame_bits bits_per_frame] = fec_encode(states, code_param, mode, payload_bits)
  ofdm_load_const;
  if code_param.data_bits_per_frame != code_param.ldpc_data_bits_per_frame
    % optionally lower the code rate by "one stuffing" - setting Nunused data bits to 1
    Nunused = code_param.ldpc_data_bits_per_frame - code_param.data_bits_per_frame;
    frame_bits = LdpcEncode([payload_bits ones(1,Nunused)], code_param.H_rows, code_param.P_matrix);
    % remove unused data bits from codeword, as they are known to the receiver and don't need to be transmitted
    frame_bits = [ frame_bits(1:code_param.data_bits_per_frame) frame_bits(code_param.ldpc_data_bits_per_frame+1:end) ];
  else
    frame_bits = LdpcEncode(payload_bits, code_param.H_rows, code_param.P_matrix);
  end
  bits_per_frame = length(frame_bits);

endfunction

function [rx_bits paritychecks] = fec_decode(states, code_param, ...
                                             payload_syms_de, payload_amps_de, ...
                                             mean_amp, EsNo)
  ofdm_load_const;
  % note ldpc_dec() handles optional lower code rate zero-stuffing
  [rx_codeword paritychecks] = ldpc_dec(code_param, mx_iter=100, demod=0, dec=0, ...
                                        payload_syms_de/mean_amp, EsNo, 
                                        payload_amps_de/mean_amp);
  rx_bits = rx_codeword(1:code_param.data_bits_per_frame);
endfunction


function [tx nclipped] = ofdm_clip(states, tx, threshold_level, plot_en=0)
  ofdm_load_const;
  tx_ = tx;
  ind = find(abs(tx) > threshold_level);
  nclipped = length(ind);
  tx(ind) = threshold_level*exp(j*angle(tx(ind)));
  if plot_en
    figure(2); clf; plot(abs(tx_(1:5*M))); hold on; plot(abs(tx(1:5*M))); hold off;
  endif
end

% two stage Hilbert clipper to improve PAPR 
function tx = ofdm_hilbert_clipper(states, tx, tx_clip_en)
  tx *= states.amp_scale;
  
  % optional compressor to improve PAPR

  nclipped = 0;
  if tx_clip_en
    if states.verbose
      printf("%f %f\n", states.clip_gain1, states.clip_gain2);
    end
    [tx nclipped] = ofdm_clip(states, tx*states.clip_gain1, states.ofdm_peak);
    
    cutoff_norm = states.txbpf_width_Hz/states.Fs;
    w_centre = mean(states.w); centre_norm = w_centre/(2*pi);
    tx = ofdm_complex_bandpass_filter(cutoff_norm, centre_norm,100,tx);
     
    % filter messes up peak levels use this to get us back to approx 16384
    tx *= states.clip_gain2;
  end

  % Hilbert Clipper 2 - remove any really low probability outliers after clipping/filtering
  % even on vanilla Tx
  [tx tmp] = ofdm_clip(states, tx, states.ofdm_peak);

  % note this is PAPR of complex signal, PAPR of real signal will be 3dB-ish larger
  peak = max(abs(tx)); RMS = sqrt(mean(abs(tx).^2));
  cpapr = 10*log10((peak.^2)/(RMS.^2));

  if states.verbose
    printf("Peak: %4.2f RMS: %5.2f CPAPR: %4.2f clipped: %5.2f%%\n",
           peak, RMS, cpapr, nclipped*100/length(tx));
  end
endfunction


% Complex bandpass filter built from low pass prototype as per src/filter.c, 
% cutoff_freq and center_freq are normalised such that cutoff_freq = 0.5 is Fs/2
function out = ofdm_complex_bandpass_filter(cutoff_freq,center_freq,n_coeffs,in)
  lowpass_coeff = fir1(n_coeffs-1, cutoff_freq);
  k = (0:n_coeffs-1);
  bandpass_coeff = lowpass_coeff .* exp(j*2*pi*center_freq*k);
  out = filter(bandpass_coeff,1,in);
endfunction


% Complex bandpass filter for Rx - just used on the very low SNR modes to help 
% with acquisition
function [rx delay_samples] = ofdm_rx_filter(states, mode, rx)
  delay_samples = 0;
  if strcmp(mode,"datac4") || strcmp(mode,"datac13")
    w_centre = mean(states.w); centre_norm = w_centre/(2*pi);
    n_coeffs = 100;
    cutoff_Hz = 400; cutoff_norm = cutoff_Hz/states.Fs;
    rx = ofdm_complex_bandpass_filter(cutoff_norm,centre_norm,n_coeffs,rx);
    delay_samples = n_coeffs/2;
  end
endfunction


% returns an unpacked CRC16 (array of 16 bits) calculated from an array of unpacked bits 
function unpacked_crc16 = crc16_unpacked(unpacked_bits)
  % pack into bytes
  mod(length(unpacked_bits),8);
  assert(mod(length(unpacked_bits),8) == 0);
  nbytes = length(unpacked_bits)/8;
  mask = 2 .^ (7:-1:0);
  for i=1:nbytes
    st = (i-1)*8 + 1; en = st+7;
    bytes(i) = sum(mask .* unpacked_bits(st:en));
  end
  crc16_hex = crc16(bytes);
  crc16_dec = [hex2dec(crc16_hex(1:2)) hex2dec(crc16_hex(3:4)) ];
  unpacked_crc16 = [];
  for b=1:length(crc16_dec)
    unpacked_crc16 = [unpacked_crc16 bitand(crc16_dec(b), mask) > 0];
  end
endfunction