#include "raw_proto.h"

/* accum_coherence_spectrum *************************************************
* INPUT:    coh     pointer to start address of complex coherence data array
*           pvrx    pointer to start of array of complex time series
*           ch      pointer to current raw data header
*           ip      input parameters from raw_main.
*           angle   last_angle array, for phase angle unwrapping
* PURPOSE:
* accumulates coherence information from pvrx in coh.  The coh array
* is complex since coherence is complex.  The coh02 array contains two sub
* arrays, coherences V0.V2* & V1.V0*
****************************************************************************/
void accum_coherence_spectrum(ComplexFloat *coh02, ComplexFloat *pvrx0,
			      HarrisHeader *ch, InputParms *ip, float *angle)
{
  ComplexFloat *last_pv, *pvrx1, *pvrx2, *coh10, temp;
  float phase,mag;
  long int rx_step;
  int i, imax;
  
  /* configured for 3 receivers only. compute coherence for all ranges */
  if((imax = ch->nchanls/2) != 3 || (ip->pm != 3 && ip->pm != 4))
  {
    printf("spectrum: accum_coherence_spectrum: invalid mode or #of rxs.  Exiting.\n",imax);
    exit(-1);
  }
  
  /* size of receiver data field. */
  rx_step = ch->tu_ng * ip->fftwdth;

  pvrx1 = pvrx0 + rx_step;
  pvrx2 = pvrx1 + rx_step;
  coh10 = coh02 + rx_step;
  
  if(ip->pm == 3)
  {
    /* compute complex phase for all ranges and doppler bins. */
    /* cmag takes conjugate of first number */
    
    for(last_pv = pvrx0 + rx_step;pvrx0 < last_pv;pvrx0++,pvrx1++,pvrx2++)
    {
      /* horizontal coherence */
      temp = cphase(pvrx2,pvrx0);
      *coh02++ = cadd(coh02,&temp);
      
      /* vertical coherence */
      temp = cphase(pvrx0,pvrx1);
      *coh10++ = cadd(coh10,&temp);
    }
    
    return;
  }
  
  if(ip->pm == 4)
  {
    /*
     * compute sums of phase angles in real part of coherence variable
     * and sums of squares of phase angles in imaginary part
     */
    /* cmag takes conjugate of first number */
    
    for(last_pv = pvrx0 + rx_step;
	pvrx0 < last_pv;
	pvrx0++,pvrx1++,pvrx2++,coh02++,coh10++)
    {
      /* horizontal phase */
      temp = cmag(pvrx2,pvrx0);
      phase = unwrap(&temp,angle++);
      coh02->r += phase;
      coh02->i += phase * phase;
      
      /* vertical phase */
      temp = cmag(pvrx0,pvrx1);
      phase = unwrap(&temp,angle++);
      coh10->r += phase;
      coh10->i += phase * phase;
    }
  }
}

/* accum_power_spectrum *****************************************************
* INPUT:    spec    pointer to start address of complex spectral data array
*           pvrx    pointer to start of array of complex time series
*           ch      pointer to current raw data header
* PURPOSE:
* accumulates power spectral information from pvrx in spec.  The spec array
* is complex to allow for double pulse phase information.
****************************************************************************/
void accum_power_spectrum(ComplexFloat *spec, ComplexFloat *pvrx, 
			  HarrisHeader *ch, InputParms *ip)
{
  ComplexFloat *last_pv,*pv_tau,temp;
  long int last_bin, tau_offset, barker_offset, sample_bins;
  int i, imax;
  
  /* only mode 3 requires accumulated power for all 3 rxs */
  imax = (ip->pm == 3) ? ch->nchanls/2: 1;
  
  /* range data, in frequency bins, that will be of interest */
  sample_bins = (ch->tu_ng - ch->n_tx_smp - ch->n_nx_smp) * ip->fftwdth;
  
  /* range data, in frequency bins, left unused at end of double pulse */
  tau_offset = ch->tu_tau / ch->tu_si * ip->fftwdth;
  
  /* range data, in frequency bins, left unused at end of barker code */
  barker_offset = (!ch->nbaudp) ? 0: (ch->nbaudp - 1) * ip->fftwdth;
  
#ifdef DEBUG_SPE
  printf("spectrum: accum_power_spectrum: sb %ld tau_os %ld barker_os %ld imax %d\n",sample_bins,tau_offset,barker_offset,imax);
#endif /* DEBUG_SPE */
  
  for(i=0;i<imax;i++)
  {
    /* first do transmitter gates if any */
    last_bin = ch->n_tx_smp * ip->fftwdth;
    for(last_pv = pvrx + last_bin;pvrx < last_pv;pvrx++)
    {
      temp = cmag(pvrx,pvrx);
      *spec++ = cadd(spec,&temp);
    }
    
    last_bin = sample_bins - tau_offset - barker_offset;
    for(last_pv = pvrx + last_bin; pvrx < last_pv; pvrx++)
    {
      temp = cmag(pvrx,pvrx + tau_offset);
      *spec++ = cadd(spec,&temp);
    }
    /* the remaining offset bins are useless, skip to noise gates */
    pvrx += tau_offset + barker_offset;
    spec += tau_offset + barker_offset;
    
    /* now do noise gates if any */
    last_bin = ch->n_nx_smp * ip->fftwdth;
    for(last_pv = pvrx + last_bin; pvrx < last_pv; pvrx++)
    {
      temp = cmag(pvrx,pvrx);
      *spec++ = cadd(spec,&temp);
    }
  }
}

/* barker ***********************************************
* INPUT:    bauds    number of bauds in barker code
*           bk       array to store barker phase code
* PURPOSE:
* generate appropriate barker code sequence or
* exit program if length not an allowed case.
********************************************************/
void barker(long int bauds, int *bk)
{
#ifdef DEBUG_SPE
  printf(
    "spectrum: barker: Finding barker sequence for %ld baud code.\n",bauds);
#endif /* DEBUG_SPE */
  
  switch(bauds)
  {
  case 3:
    *bk++ = 1;
    *bk++ = 1;        /* 1 = add, 0 = subtract */
    *bk   = 0;
    break;
  case 5:
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 0;
    *bk   = 1;
    break;
  case 7:
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 0;
    *bk++ = 0;
    *bk++ = 1;
    *bk   = 0;
    break;
  case 11:
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 0;
    *bk++ = 0;
    *bk++ = 0;
    *bk++ = 1;
    *bk++ = 0;
    *bk++ = 0;
    *bk++ = 1;
    *bk   = 0;
    break;
  case 13:
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 0;
    *bk++ = 0;
    *bk++ = 1;
    *bk++ = 1;
    *bk++ = 0;
    *bk++ = 1;
    *bk++ = 0;
    *bk   = 1;
    break;
  default:
    printf_exit("spectrum: barker: No provision made for %ld baud code.\n",
                bauds);
    break;
  }
}

/* do_barker ************************************************************
* INPUT:    pvrx    array of time sequence of complex number structs
*           ch      current Harris raw data header
*           bk      array of barker phase codings (1-add,0-sub)
* PURPOSE:
* performs n-baud barker decoding on a array of time sequences.
************************************************************************/
void do_barker(ComplexFloat *pvrx, HarrisHeader *ch, InputParms *ip, int *bk)
{
  ComplexFloat *last_pv;
  long int last_bkgate,first_bkgate,nb,offset,endskip;
  int i,imax;
  
  imax = (ip->pm == 3 || ip->pm == 4) ? ch->nchanls/2: 1;
  
#ifdef DEBUG_BARK
  printf("spectrum: do_barker: Using %ld baud decode on %d rx\n",
	 ch->nbaudp,imax);
#endif /* DEBUG_BARK */
  
  first_bkgate = ch->n_tx_smp * ip->fftwdth;
  last_bkgate = (ch->tu_ng - ch->n_nx_smp - ch->nbaudp + 1) * ip->fftwdth;
  endskip = (ch->n_nx_smp + ch->nbaudp - 1) * ip->fftwdth;
  
  if(last_bkgate < first_bkgate)
    printf_exit(
      "spectrum: do_barker: Not enough ranges for BK decoding, %ld ranges\n",
	 (last_bkgate-first_bkgate)/ip->fftwdth);
  
  for(i=0;i<imax;i++)
  {
    last_pv = pvrx + last_bkgate;
    for(pvrx += first_bkgate; pvrx < last_pv; pvrx++)
    {
      offset = ip->fftwdth;
      for(nb = 1;nb<ch->nbaudp;nb++)
      {
	*pvrx = (!bk[nb]) ?
	  csub(pvrx,pvrx + offset)
	  :
	  cadd(pvrx,pvrx + offset);
	offset += ip->fftwdth;
      }
#ifdef DEBUG_BARK
      if(!i && pvrx == last_pv - ip->fftwdth)
      {
	printf("spectrum: do_barker: Barker seq: ");
	for(nb = 0;nb<ch->nbaudp;nb++)
	{
	  printf("%d ",bk[nb]);
	}
	printf("\n");
      }
#endif /* DEBUG_BARK */
    }
    /* skip to start of next rx in case of mode 3. */
    pvrx += endskip;
  }
}

/* zero_array ***************************************
* INPUT:    arr      an array of complex numbers
*           count    how many complex numbers in arr
* PURPOSE:
* zeros an array of size count of complex numbers
****************************************************/
void zero_array(ComplexFloat *arr, long int count)
{
  ComplexFloat *p_last;
  
#ifdef THINK_C
  ComplexFloat cfzero = {0.0,0.0};
#else
  ComplexFloat cfzero;
  
  cfzero.r = cfzero.i = 0.0;
#endif /* THINK_C */
  
#ifdef DEBUG_SPE
  printf("spectrum: zero_array: Zeroing array of %ld complex numbers.\n",count);
#endif /* DEBUG_SPE */
  
  for(p_last = arr + count;arr < p_last;arr++)
    *arr = cfzero;
}

/* compute_spectra **********************************
* INPUT:    ph      processed data header           *
*           pd      float array to store spectra    *
*           spec    Complex spectral array          *
* PURPOSE:                                          *
* place spectra and header into a continous stream  *
* single pulse spectra will be reduced to floats,   *
* double pulse spectra will remain complex floats   *
* no dc averaging is done at this point.            *
****************************************************/
void compute_spectra(ProcHeader *ph, float *pd, ComplexFloat *spec)
{
  float *pspec,temp;
  long int i,j,last_r;
  ProcHeader *tmpPh;
  ComplexFloat *pnoise;
  
  /* find address where processed data header will be written */
  tmpPh = (ProcHeader *) pd;
  
  /* put start of spectral data at end of processed data header */
  pspec = (float *)((char *)pd + ph->size);
  
  /* compute noise level to store in processed data header */
  ph->bin_noise = 0.0; pnoise = &spec[ph->tu_ng - ph->n_nx_smp];
  for(i=0;i<ph->n_nx_smp;i++)
  {
    temp = compute_sp_power(pnoise,ph,1);
    if(!i)
    {
      ph->bin_noise = temp;
    }
    else
    {
      if(ph->bin_noise > temp) ph->bin_noise = temp;
    }
    pnoise += ph->fftwdth;
  }
  
  /*
   *  take the minimum of the 50th percentile noise bins
   */
  
  /*
   * store spectra for desired bins.
   * No transmitter gate or noise gate spectra will be saved.
   */
  ph->inc_tx = ph->inc_nx = 0;
  spec += (ph->n_tx_smp + ph->first_w) * ph->fftwdth;
  
  for(i=ph->first_w;i<ph->last_w;i++)
  {
    for(j=0;j<ph->fftwdth;j++)
    {
      *pspec++ = spec->r;
      /* save imaginary part only if double pulse data */
      if(ph->tu_tau) *pspec++ = spec->i;
      spec++;
    }
  }
  /* how many bytes of data and how many bytes in record */
  ph->nbytes = ph->fftwdth * 
    (ph->tu_tau ? 2:1) * 
      (ph->last_w - ph->first_w) *
	sizeof(float);
  ph->totnbyts = ph->nbytes + sizeof(ProcHeader);
  
  /* store processed data header in float stream */
  *tmpPh = *ph;
}

/* compute_coherence ********************************
* INPUT:    ph      processed data header           *
*           pd      float array to store spectra    *
*           spec0   Complex spectral array          *
*           coh02   Complex cross-coherence array   *
* PURPOSE:                                          *
* place spectral arrays and header into a continous *
* stream single pulse spectra will be reduced to    *
* floats, double pulse spectra will remain complex  *
* floats no dc averaging is done at this point.     *
* Coherence is self-normalized.                     *
****************************************************/
void compute_coherence(ProcHeader *ph, float *pd, ComplexFloat *spec0, ComplexFloat *coh02)
{
  float *pspec, *pcoh02, *pcoh10, temp;
  long int i,j,last_r,rx_step,rx_wrt,rx_first;
  ProcHeader *tmpPh;
  ComplexFloat *pnoise0, *spec1, *spec2, *coh10;
  float np0;
  
  /* find address where processed data header will be written */
  tmpPh = (ProcHeader *) pd;
  
  /* put start of spectral data at end of processed data header */
  rx_wrt = (ph->last_w - ph->first_w);
  pspec = (float *)((char *)pd + ph->size);
  pcoh02 = pspec + ph->fftwdth * (ph->tu_tau ? 2:1) * rx_wrt;
  pcoh10 = pcoh02 + ph->fftwdth * 2 * rx_wrt;
  
  rx_step = ph->tu_ng * ph->fftwdth;
  spec1 = spec0 + rx_step;
  spec2 = spec1 + rx_step;
  coh10 = coh02 + rx_step;
  
  /* compute noise level to store in processed data header */
  ph->bin_noise = 0.0;
  np0 = 0.0;
  pnoise0 = &spec0[ph->n_tx_smp];

  /*
  pnoise0 = &spec0[ph->tu_ng - ph->n_nx_smp];
  pnoise1 = &spec1[ph->tu_ng - ph->n_nx_smp];
  pnoise2 = &spec2[ph->tu_ng - ph->n_nx_smp];
  */
  for(i=0;i<ph->n_nx_smp;i++)
  {
    temp = compute_sp_power(pnoise0,ph,1);
    if(!i)
    {
      np0 = temp;
    }
    else
    {
      if(np0 > temp) np0 = temp;
    }
    pnoise0 += ph->fftwdth;
  }
  
  /*
   *  take the minimum 50th percentile noise bins
   */
  ph->bin_noise = np0;
  
  /*
   * store spectra & coherence for desired bins.
   * No transmitter gate or noise gate data will be saved.
   */
  ph->inc_tx = ph->inc_nx = 0;
  rx_first = (ph->n_tx_smp + ph->first_w);
  spec0 += rx_first * ph->fftwdth;
  spec1 += rx_first * ph->fftwdth;
  spec2 += rx_first * ph->fftwdth;
  coh02 += rx_first * ph->fftwdth;
  coh10 += rx_first * ph->fftwdth;
  
  if(ph->mode == 3)
  {
    for(i=ph->first_w;i<ph->last_w;i++)
    {
      for(j=0;j<ph->fftwdth;j++)
      {
	*pspec++ = spec0->r;
	/* save imaginary part only if double pulse data */
	if(ph->tu_tau) *pspec++ = spec0->i;
	
	temp = ph->nave;
	*pcoh02++ = coh02->r / temp;
	*pcoh02++ = coh02->i / temp;
	
	temp = ph->nave;
	*pcoh10++ = coh10->r / temp;
	*pcoh10++ = coh10->i / temp;
	
	spec0++; spec1++; spec2++; coh02++; coh10++;
      }
    }
  }
  else if(ph->mode == 4)
  {
    for(i=ph->first_w;i<ph->last_w;i++)
    {
      for(j=0;j<ph->fftwdth;j++)
      {
	*pspec++ = spec0->r;
	/* save imaginary part only if double pulse data */
	if(ph->tu_tau) *pspec++ = spec0->i;
	
	/*
	 * real part is Sum Theta
	 * imaginary part is Sum Theta Squared
	 */
	*pcoh02++ = coh02->r;
	*pcoh02++ = coh02->i;
	
	*pcoh10++ = coh10->r;
	*pcoh10++ = coh10->i;
	
	spec0++; coh02++; coh10++;
      }
    }
  }
  
  /* how many bytes of data and how many bytes in record */
  ph->nbytes = ph->fftwdth * (ph->tu_tau ? 6:5) * rx_wrt * sizeof(float);
  ph->totnbyts = ph->nbytes + sizeof(ProcHeader);
  
  /* store processed data header in float stream */
  *tmpPh = *ph;
}