testsd.c

来自「speech signal process tools」· C语言 代码 · 共 834 行 · 第 1/2 页

    case LONG:
    case SHORT:
    case DOUBLE:
      data_format = num_type;
      break;
    case DOUBLE_CPLX:
    case FLOAT_CPLX:
    case LONG_CPLX:
    case SHORT_CPLX:
    case BYTE_CPLX:
      data_format = num_type;
      c_flag = 1;
      break;
    default:
      Fprintf(stderr, "testsd: unsupported data type\n");
      exit(1);
    }
  }

  /* holdover -i option can still override all for short type */

  if (i_flag) {
    if (c_flag) 
      data_format = SHORT_CPLX;
    else 
      data_format = SHORT;
  }

  /*compute value for max_value in header*/
  switch (type) {
  case SINE:
    maxval = level; /*true for both real and complex*/
    break;
  case SQUARE:
  case TRIANGLE:
  case SAWTOOTH:
  case CONSTANT:
  case PULSES:
  case UNIFORM:
    /*sqrt(2) * level for complex*/
    maxval = (c_flag ? level * 1.414213562 : level);
    break;
  case ASCII: /*for ascii, we'll know when we get there; for 
		gauss, we just don't know*/
  case GAUSS:
    maxval = 0;
    break;
  default:
    Fprintf(stderr, "testsd: unknown test data type\n");
    exit(1);
    break;
  }

/*
 *  create and write FEA_SD header
 */
  if (debug_level) 
    Fprintf(stderr, "testsd: creating ESPS FEA_SD header\n");

  oh = new_header(FT_FEA);
  start_time = 0.0;
  (void) init_feasd_hd(oh, data_format, (int) 1, 
		       &start_time, (int) 0, (double) srate);
  (void) strcpy (oh->common.prog, program);
  (void) strcpy (oh->common.vers, version);
  (void) strcpy (oh->common.progdate, date);
  add_comment (oh, get_cmd_line(argc,argv));
  oh->common.tag = NO;

  *add_genhd_e("test_type", (short *) NULL, 1, sig_types, oh) = type;

  if (noise) *add_genhd_l("seed", (long *) NULL, 1, oh) = seed;

  if (periodic) {
    *add_genhd_f("frequency", (float *) NULL, 1, oh) = freq;
    *add_genhd_f("phase", (float *) NULL, 1, oh) = phase;
  }

  if (!a_flag)
    *add_genhd_f("level", (float *) NULL, 1, oh) = level;

  *add_genhd_d("max_value", (double *) NULL, 1, oh) = maxval;

  if (sweep_rate != 0) 
    *add_genhd_f("sweep_rate", (float *) NULL, 1, oh) = sweep_rate;

  if (decay_time != 0) 
    *add_genhd_f("decay_time", (float *) NULL, 1, oh) = decay_time;

  if (debug_level) 
    Fprintf(stderr, "testsd: writing ESPS FEA_SD header\n");

/* if we are not converting ASCII data, we write the header now; this is 
 * because the output is buffered in the case of non-ascii data.  For 
 * ascii data, the output is not buffered; memory is allocated for the
 * entire data in atoarray(); since atoarray() also returns the max value, 
 * we delay writing the header in the case of ASCII output*/

  if (type != ASCII) 
    (void) write_header(oh, ostrm);

/*
 *  allocate FEA_SD record (we use FLOAT for the in-memory records) 
 */

  if (c_flag) {
    sd_rec = allo_feasd_recs(oh, DOUBLE_CPLX, (long) BUFLEN,(char *) NULL, NO);
    if (sd_rec == NULL) {
      Fprintf(stderr, "testsd: can't allocate FEA_SD records\n");
      exit(1);
    }
    csdata = (double_cplx *) sd_rec->data;
  }
  else {
    sd_rec = allo_feasd_recs(oh, DOUBLE, (long) BUFLEN, (char *) NULL, NO);
   if (sd_rec == NULL) {
      Fprintf(stderr, "testsd: can't allocate FEA_SD records\n");
      exit(1);
    }    
    fsdata = (double *) sd_rec->data;
  }

  if (debug_level) {

    Fprintf(stderr, "testsd: allocated size %ld %s buffer\n", 
	    BUFLEN, (c_flag ? "DOUBLE_CPLX" : "DOUBLE"));
    Fprintf(stderr, 
	    "testsd: output is %d points of signal type %s\n", 
	    points, sig_types[type]);
    if (periodic)
      Fprintf(stderr, 
      "   freq = %g, phase = %g, sweep_rate = %g, decay_time = %g, level = %g\n", 
	      freq, phase, sweep_rate, decay_time, level);    
    if (type == UNIFORM) 
      Fprintf(stderr, "\tuniform noise in interval [-%g,%g]\n", 
	      level, level);
    if (type == GAUSS)
      Fprintf(stderr, "\tnoise with rms level %g\n", level);
  }
  
/*
 * generate test data
 */

  swept = (sweep_rate != 0.0);
  decay = (decay_time != 0.0);

  /* for freq = A and sweep_rate = B, the periodic functions essentially 
   * are f(TWOPI*(A + Bt/2)t), since the "instantaneous 
   * frequency" is A + Bt (derivative of phase); hence for convenience 
   * in the computations, we divide sweep_rate by 2 here; note that the 
   * variable ifreq below isn't really the instantaneous frequency, but
   * rather a sort of average over the entire waveform to that point
   */

  if (swept) 
    sweep_rate /= 2;

  points_left = points;
  block_start = -BUFLEN;

  while (points_left > 0) {
    
    c_points = ( points_left >= BUFLEN ? BUFLEN : points_left);
    points_left -= c_points;
    block_start += BUFLEN;
    if (debug_level > 1) 
      Fprintf(stderr, "testsd: generating %ld points\n", c_points);
    
    switch (type) {
    case SINE:
      /*compute time interval between SD points*/
      delta = 1 / srate;
      /*fill in data array*/
      if (c_flag) {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  csdata[i].real = ampl * cos(TWOPI * j * ifreq * delta + phase);
	  csdata[i].imag = ampl * sin(TWOPI * j * ifreq * delta + phase);
	}
      }
      else {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  value = ampl * sin(TWOPI * j * ifreq * delta + phase);
	  fsdata[i] = value;
	}
      }
      break;
    case SQUARE:
      delta = 1 / srate;
      for (i = 0; i < c_points; i++) {
	j = i + block_start;
	ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	value = ampl * square( j * ifreq * delta + TWOPI *phase);
	if (c_flag) 
	  csdata[i].real = csdata[i].imag = value;
	else
	  fsdata[i] = value;
      }
      break;
    case TRIANGLE:
      delta = 1 / srate;
      for (i = 0; i < c_points; i++)  {
	j = i + block_start;
	ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	value = ampl * triangle( j * ifreq * delta + TWOPI * phase);
	if (c_flag) 
	  csdata[i].real = csdata[i].imag = value;
	else
	  fsdata[i] = value;
      }
      break;
    case SAWTOOTH:
      delta = 1 / srate;
      for (i = 0; i < c_points; i++) {
	j = i + block_start;
	ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	value = level * sawtooth( j * ifreq * delta + TWOPI * phase);
	if (c_flag) 
	  csdata[i].real = csdata[i].imag = value;
	else
	  fsdata[i] = value;
      }
      break;
    case PULSES:
      delta = 1 / srate;
      for (i = 0; i < c_points; i++) {
	j = i + block_start;
	ifreq = (swept ? freq + sweep_rate * j * delta : freq);
	ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	cycles = j * ifreq * delta + phase/TWOPI;
	if ((long) cycles > old_cycles) {
	  value = ampl;
	  old_cycles = (long) cycles;
	}
	else 
	  value = 0;

	if (cycles == 0) /*special case for first point*/
	  value = ampl;

	if (c_flag) 
	  csdata[i].real = csdata[i].imag = value;
	else
	  fsdata[i] = value;
      }
      break;
    case GAUSS:
      /*compute time interval between SD points*/
      delta = 1 / srate;
      /*fill in data array*/
      if (c_flag) {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  rms_amp = (decay ? level * exp(- j * delta / decay_time) : level);
	  value = rms_amp * gauss();
	  csdata[i].real = value;
	  value = rms_amp * gauss();
	  csdata[i].imag = value;
	}
      }
      else {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  rms_amp = (decay ? level * exp(- j * delta / decay_time) : level);
	  value = rms_amp * gauss();
	  fsdata[i] = value;
	}
      }
      break;
    case UNIFORM:
      /*compute time interval between SD points*/
      delta = 1 / srate;
      /*fill in data array*/
      if (c_flag) {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  value = 2.0*ampl*((float)random() / BIGRAND - .5);
	  csdata[i].real = value;
	  value = 2.0*ampl*((float)random() / BIGRAND - .5);
	  csdata[i].imag = value;
	}
      }
      else {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  value = 2.0*ampl*((float)random() / BIGRAND - .5);
	  fsdata[i] = value;
	}
      }
      break;
    case CONSTANT:
      if (c_flag) {
	for (i = 0; i < c_points; i++) {
	  j = i + block_start;
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  csdata[i].real = csdata[i].imag = ampl;
	}
      }
      else {
	maxval = level;
	for (i = 0; i < c_points; i++)  {
	  ampl = (decay ? level * exp(- j * delta / decay_time) : level);
	  fsdata[i] = ampl;
	}
      }
      break;
    case ASCII:  
      if (c_flag) 
	sd_rec->data =  atoarray(instrm, DOUBLE_CPLX, &c_points, &maxval);
      else
	sd_rec->data =  atoarray(instrm, DOUBLE, &c_points, &maxval);
      
      if (debug_level) 
	Fprintf(stderr, 
	  "testsd: read %ld %s points from ASCII file %s with maxval %g\n", 
		c_points, (c_flag ? "complex" : "real"), infile, maxval);

      /* note that the allocated sd_rec->data is now irrelevant; furthermore
       * we must change num_records to be the amount returned by atoarray*/

      sd_rec->num_records = c_points; 

      /* now we can fill in max_value and write the header */
      
      *add_genhd_d("max_value", (double *) NULL, 1, oh) = maxval;

      (void) write_header(oh, ostrm);

      /*make sure loop stops since atoarray reads it all in*/


      points_left = 0;
      break;
    default:
      Fprintf(stderr, "testsd: unknown test data type\n");
      exit(1);
      break;
    }
    
    /* note that put_feasd_recs does the appropriate conversion 
     * for the type specified in the header */
    
    (void) put_feasd_recs(sd_rec, 0L,  c_points, oh, ostrm);
    
    tot_points += c_points;

    /*end of output loop*/
    
  }

/*
 * put info in ESPS common
 */

  if (strcmp(sd_file, "<stdout>") != 0) {
    (void) putsym_s("filename", sd_file);
    (void) putsym_s("prog","testsd");
    (void) putsym_i("start", (int) 1);
    (void) putsym_i("nan", (int) tot_points);
  }

/*
 * clean up and exit
 */
  exit(0);
}

double
square (angle)
double angle;
{
  angle = fmod(angle, 1.0); 

  if (angle < 0.5)
    return 1.0;

  return -1.0;
}

double
triangle(angle)
double angle;
{
  angle = fmod(angle, 1.0); 
 
  if (angle <= 0.25)
    return (4 * angle);

  if ((angle > 0.25) && (angle <= 0.75))
    return (2 - 4 * angle);

  return (-4.0 + 4 * angle);
}

double
sawtooth(angle)
double angle;
{
  return(fmod(angle, 1.0)); 
}