/*-------------- Telecommunications & Signal Processing Lab ---------------
                             McGill University

Routine:
  void RSexpTime (double t, double Ns, long int *n, double *ds, double *dsr)

Purpose:
  Express a value in terms of residues

Description:
  This procedure expresses a given value t in terms of residues n, ds and dsr,
    t = n + (dsr + ds) / Ns,
  where n is an integer value, ds is a positive value less than Ns, and dsr
  is a positive fractional value, 0 <= dsr < 1.  If Ns is an integer value,
  then ds will also be integer.  If Ns is not integer, then dsr is set to
  zero and ds will have a fractional component.

Parameters:
   -> double t
      Input value to be represented
   -> double Ns
      Input value used to perform the modulo operations
  <-  long int *n
      Integral part of the representation
  <-  double *ds
      Output value used to represent the input value, 0 <= ds < Ns.
  <-  double *dsr
      Output value used to represent the input value, 0 <= dsr < 1.

Author / revision:
  P. Kabal  Copyright (C) 1996
  $Revision: 1.3 $  $Date: 1996/07/16 15:53:22 $

-------------------------------------------------------------------------*/

static char rcsid[] = "$Id: RSexpTime.c 1.3 1996/07/16 AFsp-V2R1 $";

#include <math.h>
#include <libtsp.h>
#include "ResampAudio.h"

#define EPS	1E-8

/*
   The value t points to the sample position for an output value.  For the
   default filter design, the filter delay compensates (to an integral multiple
   of 1/Ir) for the fractional part of the sample offset.  The value Ns is
   set to be a multiple of Ir if possible.  In that case, the fractional sample
   offset can be expressed as
     frac(t) = (dsr + ds) / Ns ,
   where ds is integer and dsr is zero.  Due to numerical roundoff in
   calculating t, dsr can be eps or 1-eps, where eps is a small number.

   Here we "force" dsr to be zero if it is near zero or one, compensating ds
   and n if necessary for the change in dsr.
*/

void
RSexpTime (t, Ns, n, ds, dsr)

     double t;
     double Ns;
     long int *n;
     double *ds;
     double *dsr;

{
  long int nv;
  double dsv;
  double dsrv;
  
  nv = (long int) floor (t);		/* floor, long int */

  /* Integer Ns */
  if (Ns == floor (Ns)) {
    dsv = floor (Ns * (t - nv));	/* double, integer-valued */
    dsrv = Ns * (t - nv) -  dsv;	/* double */

    /* Move dsr if it is nearly an integer */
    if (dsrv < EPS)
      dsrv = 0.0;
    else if (dsrv > 1.0-EPS) {
      dsrv = 0.0;
      dsv = dsv + 1.0;
      if  (dsv >= Ns) {
	dsv = 0.0;
	++nv;
      }
    }
  }

  /* Non-integer Ns */
  else {
    dsv = Ns * (t - nv);
    dsrv = 0.0;
  }

  *n = nv;
  *ds = dsv;
  *dsr = dsrv;

  return;
}