/*************************************************************************/
/*                                                                       */
/*                Centre for Speech Technology Research                  */
/*                     University of Edinburgh, UK                       */
/*                         Copyright (c) 1997                            */
/*                        All Rights Reserved.                           */
/*                                                                       */
/*  Permission is hereby granted, free of charge, to use and distribute  */
/*  this software and its documentation without restriction, including   */
/*  without limitation the rights to use, copy, modify, merge, publish,  */
/*  distribute, sublicense, and/or sell copies of this work, and to      */
/*  permit persons to whom this work is furnished to do so, subject to   */
/*  the following conditions:                                            */
/*   1. The code must retain the above copyright notice, this list of    */
/*      conditions and the following disclaimer.                         */
/*   2. Any modifications must be clearly marked as such.                */
/*   3. Original authors' names are not deleted.                         */
/*   4. The authors' names are not used to endorse or promote products   */
/*      derived from this software without specific prior written        */
/*      permission.                                                      */
/*                                                                       */
/*  THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK        */
/*  DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING      */
/*  ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT   */
/*  SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE     */
/*  FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES    */
/*  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN   */
/*  AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,          */
/*  ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF       */
/*  THIS SOFTWARE.                                                       */
/*                                                                       */
/*************************************************************************/
/*                      Author :  Paul Bagshaw                           */
/*                      Date   :  1993                                   */
/*************************************************************************/
/*                                                                       */
/* The above copyright was given by Paul Bagshaw, he retains             */
/* his original rights                                                   */
/*                                                                       */
/*************************************************************************/
 /****************************************************************************
 *                                                                           *
 * Pitch Determination Algorithm.                                            *
 *                                                                           *
 * Super Resolution Pitch Determinator with No Headers (SRPD_HD).            *
 *                                                                           *
 * Analysis synchronised with cepstral analysis, pitch biasing option, and   *
 * optimised for minimum gross pitch errors and accurate voiced/unvoiced     *
 * classification. All known bugs resolved!                                  *
 *                                                                           *
 * 4th February 1992:                                                        *
 * Additional option [-w] added to give an artificial frame length, thus     *
 * allowing the output data to be synchronised with other signal processing  *
 * algorithms such as cepstral analysis and formant tracking.                *
 *                                                                           *
 * Y. Medan, E. Yair, and D. Chazan, "Super resolution pitch determination   *
 * of speech signals," IEEE Trans. Signal Processing Vol.39 No.1             *
 * pp.40-48 (1991).                                                          *
 *                                                                           *
 * Implementation by Paul Bagshaw, Centre for Speech Technology Research,    *
 * University of Edinburgh, 80 South Bridge, Edinburgh EH1 1HN.              *
 *                                                                           *
 *****************************************************************************/

/************************
 * include header files *
 ************************/

#include <cmath>
#include <cstdlib>
#include <iostream>
#include "srpd.h"
#include "EST_cutils.h"
#include "EST_Wave.h"

#ifndef MAXSHORT
#define MAXSHORT 32767
#endif

void super_resolution_pda (struct Srpd_Op *paras, SEGMENT_ seg,
			   CROSS_CORR_ *p_cc, STATUS_ *p_status)
{

  static int zx_lft_N, zx_rht_N;
  static double prev_pf = BREAK_NUMBER;

  int n, j, k, N0 = 0, N1, N2, N_, q, lower_found = 0, score = 1, apply_bias;
  int x_index, y_index, z_index;
  int zx_rate = 0, zx_at_N0 = 0, prev_sign;
  int seg1_zxs = 0, seg2_zxs = 0, total_zxs;
  short prev_seg1, prev_seg2;
  short x_max = -MAXSHORT, x_min = MAXSHORT;
  short y_max = -MAXSHORT, y_min = MAXSHORT;
  double xx = 0.0, yy = 0.0, zz = 0.0, xy = 0.0, yz = 0.0, xz = 0.0;
  double max_cc = 0.0, coefficient, coeff_weight;
  double xx_N, yy_N, xy_N, y1y1_N, xy1_N, yy1_N, beta;
  LIST_ *sig_pks_hd, *sig_pks_tl, *sig_peak, *head, *tail;
  
  sig_pks_hd = head = NULL;
  sig_pks_tl = tail = NULL;
  /* set correlation coefficient threshold */
  if (p_status->v_uv == UNVOICED || p_status->v_uv == SILENT)
    p_status->threshold = paras->Thigh;
  else /* p_status->v_uv == VOICED */
    p_status->threshold = (paras->Tmin > paras->Tmax_ratio *
        p_status->cc_max) ? paras->Tmin : paras->Tmax_ratio *
	p_status->cc_max;
  /* determine if a bias should be applied */
  if (paras->peak_tracking && prev_pf != BREAK_NUMBER &&
      p_status->v_uv == VOICED && p_status->s_h != HOLD &&
      p_status->pitch_freq < 1.75 * prev_pf &&
      p_status->pitch_freq > 0.625 * prev_pf)
    apply_bias = 1;
  else
    apply_bias = 0;
  /* consider first two segments of period n = Nmin */
  prev_seg1 = seg.data[paras->Nmax - paras->Nmin] < 0 ? -1 : 1;
  prev_seg2 = seg.data[paras->Nmax] < 0 ? -1 : 1;
  for (j = 0; j < paras->Nmin; j += paras->L) {
    /* find max and min amplitudes in x and y segments */
    x_index = paras->Nmax - paras->Nmin + j;
    y_index = paras->Nmax + j;
    if (seg.data[x_index] > x_max) x_max = seg.data[x_index];
    if (seg.data[x_index] < x_min) x_min = seg.data[x_index];
    if (seg.data[y_index] > y_max) y_max = seg.data[y_index];
    if (seg.data[y_index] < y_min) y_min = seg.data[y_index];
    /* does new sample in x or y segment represent an input zero-crossing */
    if (seg.data[x_index] * prev_seg1 < 0) {
      prev_seg1 *= -1;
      seg1_zxs++;
    }
    if (seg.data[y_index] * prev_seg2 < 0) {
      prev_seg2 *= -1;
      seg2_zxs++;
    }
    /* calculate parts for first correlation coefficient */
    xx += (double) seg.data[x_index] * seg.data[x_index];
    yy += (double) seg.data[y_index] * seg.data[y_index];
    xy += (double) seg.data[x_index] * seg.data[y_index];
  }
  /* low amplitude segment represents silence */
  if (abs (x_max) + abs (x_min) < 2 * paras->Tsilent || 
      abs (y_max) + abs (y_min) < 2 * paras->Tsilent) {
    for (q = 0; q < p_cc->size; p_cc->coeff[q++] = 0.0);
    prev_pf = p_status->pitch_freq;
    p_status->pitch_freq = BREAK_NUMBER;
    p_status->v_uv = SILENT;
    p_status->s_h = SEND;
    p_status->cc_max = 0.0;
    return;
  }
  /* determine first correlation coefficients, for period n = Nmin */
  p_cc->coeff[0] = p_status->cc_max = xy / sqrt (xx) / sqrt (yy);
  for (q = 1; q < p_cc->size && q < paras->L; p_cc->coeff[q++] = 0.0);
  total_zxs = seg1_zxs + seg2_zxs;
  prev_sign = p_cc->coeff[0] < 0.0 ? -1 : 1;
  prev_seg1 = seg.data[paras->Nmax - paras->Nmin] < 0 ? -1 : 1;
  /* iteratively determine correlation coefficient for next possible period */
  for (n = paras->Nmin + paras->L; n <= paras->Nmax; n += paras->L,
       j += paras->L) {
    x_index = paras->Nmax - n;
    y_index = paras->Nmax + j;
    /* does new samples in x or y segment represent an input zero-crossing */
    if (seg.data[x_index] * prev_seg1 < 0) {
      prev_seg1 *= -1;
      total_zxs++;
    }
    if (seg.data[y_index] * prev_seg2 < 0) {
      prev_seg2 *= -1;
      total_zxs++;
    }
    /* determine next coefficient */
    xx += (double) seg.data[x_index] * seg.data[x_index];
    yy += (double) seg.data[y_index] * seg.data[y_index];
    for (k = 0, xy = 0.0; k < n; k += paras->L)
      xy += (double) seg.data[paras->Nmax - n + k] * seg.data[paras->Nmax + k];
    p_cc->coeff[n - paras->Nmin] = xy / sqrt (xx) / sqrt (yy);
    if (p_cc->coeff[n - paras->Nmin] > p_status->cc_max)
      p_status->cc_max = p_cc->coeff[n - paras->Nmin];
    /* set unknown coefficients to zero */
    for (q = n - paras->Nmin + 1;
	 q < p_cc->size && q < n - paras->Nmin + paras->L;
	 p_cc->coeff[q++] = 0.0);
    /* is there a slope with positive gradient in the coefficients track yet */
    if (p_cc->coeff[n - paras->Nmin] > p_cc->coeff[n - paras->Nmin - paras->L])
      lower_found = 1;
    /* has new coefficient resulted in a zero-crossing */
    if (p_cc->coeff[n - paras->Nmin] * prev_sign < 0.0) {
      prev_sign *= -1;
      zx_rate++;
    }
    /* does the new coefficient represent a pitch period candidate */
    if (N0 != 0 && zx_rate > zx_at_N0) {
      add_to_list (&sig_pks_hd, &sig_pks_tl, N0, 1);
      N0 = 0;
      max_cc = 0.0;
    }
    if (apply_bias && n > zx_lft_N && n < zx_rht_N)
      coeff_weight = 2.0;
    else
      coeff_weight = 1.0;
    if (p_cc->coeff[n - paras->Nmin] > max_cc && total_zxs > 3 && lower_found) {
      max_cc = p_cc->coeff[n - paras->Nmin];
      if (max_cc * coeff_weight >= p_status->threshold) {
	zx_at_N0 = zx_rate;
	N0 = n;
      }
    }
  }
  /* unvoiced if no significant peak found in coefficients track */
  if (sig_pks_hd == NULL) {
    prev_pf = p_status->pitch_freq;
    p_status->pitch_freq = BREAK_NUMBER;
    p_status->v_uv = UNVOICED;
    p_status->s_h = SEND;
    return;
  }
  /* find which significant peak in list corresponds to true pitch period */
  sig_peak = sig_pks_hd;
  while (sig_peak != NULL) {
    yy = zz = yz = 0.0;
    for (j = 0; j < sig_peak->N0; j++) {
      y_index = paras->Nmax + j;
      z_index = paras->Nmax + sig_peak->N0 + j;
      yy += (double) seg.data[y_index] * seg.data[y_index];
      zz += (double) seg.data[z_index] * seg.data[z_index];
      yz += (double) seg.data[y_index] * seg.data[z_index];
    }
    if (yy == 0.0 || zz == 0.0)
      coefficient = 0.0;
    else
      coefficient = yz / sqrt (yy) / sqrt (zz);
    if (apply_bias && sig_peak->N0 > zx_lft_N && sig_peak->N0 < zx_rht_N)
      coeff_weight = 2.0;
    else
      coeff_weight = 1.0;
    if (coefficient * coeff_weight >= p_status->threshold) {
      sig_peak->score = 2;
      if (head == NULL) {
	head = sig_peak;
	score = 2;
      }
      tail = sig_peak;
    }
    sig_peak = sig_peak->next_item;
  }
  if (head == NULL) head = sig_pks_hd;
  if (tail == NULL) tail = sig_pks_tl;
  N0 = head->N0;
  if (tail != head) {
    xx = 0.0;
    for (j = 0; j < tail->N0; j++)
      xx += (double) seg.data[paras->Nmax - tail->N0 + j] *
	  seg.data[paras->Nmax - tail->N0 + j];
    sig_peak = head;
    while (sig_peak != NULL) {
      if (sig_peak->score == score) {
	xz = zz = 0.0;
	for (j = 0; j < tail->N0; j++) {
	  z_index = paras->Nmax + sig_peak->N0 + j;
	  xz += (double) seg.data[paras->Nmax - tail->N0 + j] *
	      seg.data[z_index];
	  zz += (double) seg.data[z_index] * seg.data[z_index];
	}
	coefficient = xz / sqrt (xx) / sqrt (zz);
	if (sig_peak == head)
	  max_cc = coefficient;
	else if (coefficient * paras->Tdh > max_cc) {
	  N0 = sig_peak->N0;
	  max_cc = coefficient;
	}
      }
      sig_peak = sig_peak->next_item;
    }
  }
  p_status->cc_max = p_cc->coeff[N0 - paras->Nmin];
  /* voiced segment period now found */
  if ((tail == head && score == 1 && p_status->v_uv != VOICED) ||
      p_cc->coeff[N0 - paras->Nmin] < p_status->threshold)
    p_status->s_h = HOLD;
  else
    p_status->s_h = SEND;
  /* find left and right boundaries of peak in coefficients track */
  zx_lft_N = zx_rht_N = 0;
  for (q = N0; q >= paras->Nmin; q -= paras->L)
    if (p_cc->coeff[q - paras->Nmin] < 0.0) {
      zx_lft_N = q;
      break;
    }
  for (q = N0; q <= paras->Nmax; q += paras->L)
    if (p_cc->coeff[q - paras->Nmin] < 0.0) {
      zx_rht_N = q;
      break;
    }
  /* define small region around peak */
  if (N0 - paras->L < paras->Nmin) {
    N1 = N0;
    N2 = N0 + 2 * paras->L;
  }
  else if (N0 + paras->L > paras->Nmax) {
    N1 = N0 - 2 * paras->L;
    N2 = N0;
  }
  else {
    N1 = N0 - paras->L;
    N2 = N0 + paras->L;
  }
  /* compensate for decimation factor L */
  if (paras->L != 1) {
    xx = yy = xy = 0.0;
    for (j = 0; j < N1; j++) {
      x_index = paras->Nmax - N1 + j;
      y_index = paras->Nmax + j;
      xx += (double) seg.data[x_index] * seg.data[x_index];
      xy += (double) seg.data[x_index] * seg.data[y_index];
      yy += (double) seg.data[y_index] * seg.data[y_index];
    }
    p_cc->coeff[N1 - paras->Nmin] = p_status->cc_max =
        xy / sqrt (xx) / sqrt (yy);
    N0 = N1;
    for (n = N1 + 1; n <= N2; n++, j++) {
      xx += (double) seg.data[paras->Nmax - n] * seg.data[paras->Nmax - n];
      yy += (double) seg.data[paras->Nmax + j] * seg.data[paras->Nmax + j];
      for (k = 0, xy = 0.0; k < n; k++)
	xy += (double) seg.data[paras->Nmax - n + k] * seg.data[paras->Nmax + k];
      p_cc->coeff[n - paras->Nmin] = xy / sqrt (xx) / sqrt (yy);
      if (p_cc->coeff[n - paras->Nmin] > p_status->cc_max) {
	p_status->cc_max = p_cc->coeff[n - paras->Nmin];
	N0 = n;
      }
    }
  }
  /* compensate for finite resolution in estimating pitch */
  if (N0 - 1 < paras->Nmin || N0 == N1) N_ = N0;
  else if (N0 + 1 > paras->Nmax || N0 == N2) N_ = N0 - 1;
  else if (p_cc->coeff[N0 - paras->Nmin] - p_cc->coeff[N0 - paras->Nmin - 1] <
	   p_cc->coeff[N0 - paras->Nmin] - p_cc->coeff[N0 - paras->Nmin + 1])
    N_ = N0 - 1;
  else
    N_ = N0;
  xx_N = yy_N = xy_N = y1y1_N = xy1_N = yy1_N = 0.0;
  for (j = 0; j < N_; j++) {
    x_index = paras->Nmax - N_ + j;
    y_index = paras->Nmax + j;
    xx_N += (double) seg.data[x_index] * seg.data[x_index];
    yy_N += (double) seg.data[y_index] * seg.data[y_index];
    xy_N += (double) seg.data[x_index] * seg.data[y_index];
    y1y1_N += (double) seg.data[y_index + 1] * seg.data[y_index + 1];
    xy1_N += (double) seg.data[x_index] * seg.data[y_index + 1];
    yy1_N += (double) seg.data[y_index] * seg.data[y_index + 1];
  }
  beta = (xy1_N * yy_N - xy_N * yy1_N) /
      (xy1_N * (yy_N - yy1_N) + xy_N * (y1y1_N - yy1_N));
  if (beta < 0.0) {
    N_--;
    beta = 0.0;
  }
  else if (beta >= 1.0) {
    N_++;
    beta = 0.0;
  }
  else
    p_status->cc_max = ((1.0 - beta) * xy_N + beta * xy1_N) /
      sqrt (xx_N * ((1.0 - beta) * (1.0 - beta) * yy_N +
		    2.0 * beta * (1.0 - beta) * yy1_N +
		    beta * beta * y1y1_N));
  prev_pf = p_status->pitch_freq;
  p_status->pitch_freq = (double) (paras->sample_freq) / (double) (N_ + beta);
  p_status->v_uv = VOICED;
  free_list (&sig_pks_hd);
  return;

}

/************* * LEVEL TWO * ************/

void add_to_list (LIST_ **p_list_hd, LIST_ **p_list_tl, int N_val, 
		  int score_val)
{

  LIST_ *new_node, *last_node;

  new_node = walloc(LIST_ ,1);
  last_node = *p_list_tl;
  new_node->N0 = N_val;
  new_node->score = score_val;
  new_node->next_item = NULL;
  if (*p_list_hd == NULL)
    *p_list_hd = new_node;
  else
    last_node->next_item = new_node;
  *p_list_tl = new_node;

}

/********************
 * define functions *
 ********************/

/************* * LEVEL ONE * ************/

void error (error_flags err_type)
{

  char prog[15]; /* program file name */

  strcpy (prog, "srpd");
  fprintf (stderr, "%s: ", prog);
  switch (err_type) {
  case CANT_WRITE:
    fprintf (stderr, "cannot write to output file");
    break;
  case DECI_FCTR:
    fprintf (stderr, "decimation factor not set");
    break;
  case INSUF_MEM:
    fprintf (stderr, "insufficient memory available");
    break;
  case FILE_ERR:
    perror ("");
    break;
  case FILE_SEEK:
    fprintf (stderr, "improper fseek () to reposition a stream");
    break;
  case LEN_OOR:
    fprintf (stderr, "artificial frame length set out of range");
    break;
  case MAX_FREQ:
    fprintf (stderr, "maximum pitch frequency value (Hz) not set");
    break;
  case MIN_FREQ:
    fprintf (stderr, "minimum pitch frequency value (Hz) not set");
    break;
  case MISUSE:
    fprintf (stderr, "usage: %s -i lpf_sample_file ", prog);
    fprintf (stderr, "-o pitch_file [options]\n");
    fprintf (stderr, "\nOptions {with default values}\n");
    fprintf (stderr, "-a form pitch_file in ascii format\n");
    fprintf (stderr, "-l 'lower pitch frequency limit' {%f (Hz)}\n",
	     DEFAULT_MIN_PITCH);
    fprintf (stderr, "-u 'upper pitch frequency limit' {%f (Hz)}\n",
	     DEFAULT_MAX_PITCH);
    fprintf (stderr, "-d 'decimation factor' {%d (samples)}\n",
	     DEFAULT_DECIMATION);
    fprintf (stderr, "-n 'noise floor (abs. amplitude)' {%d}\n",
	     DEFAULT_TSILENT);
    fprintf (stderr, "-h 'unvoiced to voiced coeff threshold' {%f}\n",
	     DEFAULT_THIGH);
    fprintf (stderr, "-m 'min. voiced to unvoiced coeff threshold' {%f}\n",
	     DEFAULT_TMIN);
    fprintf (stderr, "-r 'voiced to unvoiced coeff threshold ratio' {%f}\n",
	     DEFAULT_TMAX_RATIO);
    fprintf (stderr, "-t 'anti pitch doubling/halving threshold' {%f}\n",
	     DEFAULT_TDH);
    fprintf (stderr, "-p perform peak tracking\n");
    fprintf (stderr, "-f 'sampling frequency' {%d (Hz)}\n", DEFAULT_SF);
    fprintf (stderr, "-s 'frame shift' {%f (ms)}\n", DEFAULT_SHIFT);
    fprintf (stderr, "-w 'artificial frame length' {%f (ms)}\n",
	     DEFAULT_LENGTH);
    break;
  case NOISE_FLOOR:
    fprintf (stderr, "noise floor set below minimum amplitude");
    break;
  case SAMPLE_FREQ:
    fprintf (stderr, "attempt to set sampling frequency negative");
    break;
  case SFT_OOR:
    fprintf (stderr, "frame shift set out of range");
    break;
  case THR_DH:
    fprintf (stderr, "anti pitch doubling/halving threshold not set");
    break;
  case THR_HIGH:
    fprintf (stderr, "unvoiced to voiced coeff threshold not set");
    break;
  case THR_MAX_RTO:
    fprintf (stderr, "voiced to unvoiced coeff threshold ratio not set");
    break;
  case THR_MIN:
    fprintf (stderr, "minimum voiced to unvoiced coeff threshold not set");
    break;
  default:
    fprintf (stderr, "undefined error, %u occurred", err_type);
    break;
  }
  fprintf (stderr, "\n");
  exit (-1);

}

void initialise_parameters (struct Srpd_Op *p_par)
{
  p_par->L = DEFAULT_DECIMATION;
  p_par->min_pitch = DEFAULT_MIN_PITCH;
  p_par->max_pitch = DEFAULT_MAX_PITCH;
  p_par->shift = DEFAULT_SHIFT;
  p_par->length = DEFAULT_LENGTH;
  p_par->Tsilent = DEFAULT_TSILENT;
  p_par->Tmin = DEFAULT_TMIN;
  p_par->Tmax_ratio = DEFAULT_TMAX_RATIO;
  p_par->Thigh = DEFAULT_THIGH;
  p_par->Tdh = DEFAULT_TDH;
  p_par->make_ascii = 0;
  p_par->peak_tracking = 0;
  p_par->sample_freq = DEFAULT_SF;
  /* p_par->Nmax and p_par->Nmin cannot be initialised */
  return;

}

void initialise_structures (struct Srpd_Op *p_par, SEGMENT_ *p_seg,
     CROSS_CORR_ *p_cc)
{
  p_par->Nmax = (int) ceil((float)p_par->sample_freq / p_par->min_pitch);
  p_par->Nmin = (int) floor((float)p_par->sample_freq / p_par->max_pitch);
  p_par->min_pitch = (float)p_par->sample_freq / (float)p_par->Nmax;
  p_par->max_pitch = (float)p_par->sample_freq / (float)p_par->Nmin;

  p_seg->size = 3 * p_par->Nmax;
  p_seg->shift = (int) rint( p_par->shift / 1000.0 * (float)p_par->sample_freq );
  p_seg->length = (int) rint( p_par->length / 1000.0 * (float)p_par->sample_freq );
  p_seg->data = walloc(short,p_seg->size);

  p_cc->size = p_par->Nmax - p_par->Nmin + 1;
  p_cc->coeff = walloc(double,p_cc->size);

  return;
}


void initialise_status (struct Srpd_Op *paras, STATUS_ *p_status)
{

  p_status->pitch_freq = BREAK_NUMBER;
  p_status->v_uv = SILENT;
  p_status->s_h = SEND; /* SENT */
  p_status->cc_max = 0.0;
  p_status->threshold = paras->Thigh;
  return;

}

void end_structure_use(SEGMENT_ *p_seg, CROSS_CORR_ *p_cc)
{

  wfree (p_seg->data);
  wfree (p_cc->coeff);
  return;

}

#define BEGINNING 1
#define MIDDLE_   2
#define END       3

int read_next_segment (FILE *voxfile, struct Srpd_Op *paras, SEGMENT_ *p_seg)
{

  static int status = BEGINNING, padding= -1, tracklen = 0;

  int samples_read = 0;
  long init_file_position, offset;

  if (status == BEGINNING) {
    if (padding == -1) {
      if (fseek (voxfile, 0L, 2)) error (FILE_SEEK);
      tracklen = ((ftell (voxfile) / sizeof (short)) - p_seg->length) /
				 p_seg->shift + 1;
	cout << "track len " << tracklen;
      rewind (voxfile);
      if (paras->Nmax < p_seg->length / 2) {
	offset = (long) (p_seg->length / 2 - paras->Nmax) * sizeof (short);
	if (fseek (voxfile, offset, 1)) error (FILE_SEEK);
	padding = 0;
      }
      else {
	if ((paras->Nmax - p_seg->length / 2) % p_seg->shift != 0) {
	  offset = (long) (p_seg->shift - ((paras->Nmax - p_seg->length / 2) %
					   p_seg->shift)) * sizeof (short);
	  if (fseek (voxfile, offset, 1)) error (FILE_SEEK);
	}
	padding = (paras->Nmax - p_seg->length / 2) / p_seg->shift +
	  ((paras->Nmax - p_seg->length / 2) % p_seg->shift == 0 ? 0 : 1);
   }
  }
    cout << "padding " << padding << endl;
    if (padding-- == 0)
      status = MIDDLE_;
    else if (tracklen-- <= 0)
      return (0);
    else
      return (2);
  }
  cout << "tl  " << tracklen << endl;
  if (status == MIDDLE_) {
    if (tracklen > 0) {
      init_file_position = ftell (voxfile);
      offset = (long) (p_seg->shift * sizeof (short));
      samples_read = fread ((short *) p_seg->data, sizeof (short),
			    p_seg->size, voxfile);
      if (samples_read == p_seg->size) {
	if (fseek (voxfile, init_file_position + offset, 0)) error (FILE_SEEK);
	tracklen--;
	return (1);
      }
      else {
	status = END;
      }
    }
    else
      return (0);
  }
  if (status == END) {
    if (tracklen-- > 0)
      return (2);
    else
      return (0);
  }
  return (0);

}

int read_next_wave_segment(EST_Wave &sig, Srpd_Op *paras, SEGMENT_ *p_seg)
{
    static int status = BEGINNING, padding = -1, tracklen = 0;
    int i;
    long offset;
    static int wave_pos;
    

    //printf("read:  size %d shift %d length %d\n", p_seg->size, p_seg->shift, p_seg->length);

    if (status == BEGINNING) 
    {
	if (padding == -1) 
	{
	    tracklen = (sig.num_samples() - p_seg->length) 
		/ p_seg->shift + 1;
	    if (paras->Nmax < p_seg->length / 2) 
	    {
		offset = p_seg->length / 2 - paras->Nmax;
		wave_pos = offset;
		padding = 0;
	    }
	    else 
	    {
		if ((paras->Nmax - p_seg->length / 2) % p_seg->shift != 0) {
		    offset = p_seg->shift - ((paras->Nmax - p_seg->length / 2)%
						     p_seg->shift);
		    wave_pos = offset;
		}
		padding = (paras->Nmax - p_seg->length / 2) / p_seg->shift +
		    ((paras->Nmax - p_seg->length / 2) 
		     % p_seg->shift == 0 ? 0 : 1);
	    }
	}
	if (padding-- == 0)
	    status = MIDDLE_;
	else if (tracklen-- <= 0) {
		status = BEGINNING;
		padding = -1;
		tracklen = 0;
	    return (0);
	}
	else
	    return (2);
    }
    if (status == MIDDLE_)
    {
	if (tracklen > 0) 
	{
	    offset = p_seg->shift;
	    for (i = 0; (i < p_seg->size) && (i+wave_pos)<sig.num_samples();
		 ++i)
		p_seg->data[i] = sig.a(i + wave_pos,0);
	    for ( ; i < p_seg->size; ++i)
		p_seg->data[i] = 0;

	    if (wave_pos <= sig.num_samples())
	    {
		wave_pos += offset;
		tracklen--;
		return (1);
	    }
	    else
		status = END;
	}
	else {
		status = BEGINNING;
		padding = -1;
		tracklen = 0;
	    return (0);
	    }
    }
    if (status == END) 
    {
	if (tracklen-- > 0)
	    return (2);
	else	{
		status = BEGINNING;
		padding = -1;
		tracklen = 0;
	    return (0);
	    }
    }
    status = BEGINNING;
    padding = -1;
    tracklen = 0;
    return (0);
}

void write_track(STATUS_ status, struct Srpd_Op paras, FILE *outfile)
{
  if (paras.make_ascii) 
  {
      if (fprintf(outfile,"%7g\n",status.pitch_freq) != 8)
	  error(CANT_WRITE);
  }
  else
    if (!fwrite ((double *) &status.pitch_freq, sizeof (double), 1, outfile))
      error (CANT_WRITE);
  return;

}

void free_list (LIST_ **p_list_hd)
{

  LIST_ *next;

  while (*p_list_hd != NULL) {
    next = (*p_list_hd)->next_item;
    wfree (*p_list_hd);
    *p_list_hd = next;
  }

}