sgram.c

speech signal process tools

	    start = getsym_i("start");
	    break;
	default:
	    ERROR_EXIT("Parameter \"start\" not of integer type.")
	    break;
	}

	switch (symtype("nan"))
	{
	case ST_UNDEF:
	    nan = 0;
	    break;
	case ST_INT:
	    nan = getsym_i("nan");
	    break;
	default:
	    ERROR_EXIT("Parameter \"nan\" not of integer type.")
	    break;
	}

	end = (nan == 0) ? LONG_MAX : start + nan - 1;
	data_start_time = file_start_time + (start - 1) / sf;
    }

    if (start < 1) {  /* In case of precision problems in start_time specs... */
         start = 1;
	 if (debug_level){
	   Fprintf(stderr,"Can't start before beginning of file.\n");
	   Fprintf(stderr, "Resetting starting sample to 1.\n");
	 }
    }


    if (debug_level) Fprintf(stderr,
	"start = %ld, nan = %ld, end = %ld, data_start_time = %g\n",
	start, nan, end, data_start_time);

    /* -w window_len */

    if (!(wflag || p_mflag))
	switch (symtype("window_len"))
	{
	case ST_UNDEF:
	    window_len = def_window_len;
	    break;
	case ST_FLOAT:
	    window_len = getsym_d("window_len");
	    break;
	default:
	    ERROR_EXIT("Parameter \"window_len\" not of float type.")
	    break;
	}

    window_size = frame_len = LROUND(window_len * sf / 1000.0);

    if((win == WT_SINC) || (win == WT_SINC_C4))
      frame_len = tr_length;

    if (debug_level)
	Fprintf(stderr, "window_len = %g, frame_len = %ld\n",
		window_len, frame_len);

    /* -E pre_emphasis */

    if (!(Eflag || p_mflag))
	switch (symtype("pre_emphasis"))
	{
	case ST_UNDEF:
	    preemph = def_preemph;
	    break;
	case ST_FLOAT:
	    preemph = getsym_d("pre_emphasis");
	    break;
	default:
	    ERROR_EXIT("Parameter \"pre_emphasis\" not of float type.")
	    break;
	}

    if (debug_level)
	Fprintf(stderr, "preemph = %g\n", preemph);

    if (!(Tflag))
	switch (symtype("desired_frames"))
	{
	case ST_UNDEF:
	    break;
	case ST_INT:
	    des_frames = getsym_i("desired_frames");
	    if (des_frames > 0)
	      Tflag++;
	    break;
	default:
	    ERROR_EXIT("Parameter \"desired_frames\" not of int type.")
	    break;
	}


    /* -S step_size */

    if (!(Sflag || p_mflag || Tflag))
	switch (symtype("step_size"))
	{
	case ST_UNDEF:
	    step_size = def_step_size;
	    break;
	case ST_FLOAT:
	    step_size = getsym_d("step_size");
	    break;
	default:
	    ERROR_EXIT("Parameter \"step_size\" not of float type.")
	    break;
	}

    step = LROUND(step_size * sf / 1000.0);


    /* If -T (TI) option is used, the number of frames (number of ffts) 
     * is forced.  This facilitates looking at a narrow region with 
     * higher resolution.   Here we compute a step size that gives 
     * the (approximate) right number of frames
     */

    if (Tflag) {
	
	if (nan == 0) {
	  Fprintf(stderr, 
	      "%s: can't use -T without specifying a range start\n", ProgName);
	  Fprintf(stderr,"\tand stop with -r or -p\n");
	  exit(1);
	}

	if (Sflag && (debug_level || !zflag))
	  Fprintf(stderr, "%s: WARNING -S option overriden by -T\n", ProgName);

	/* computed fixed step size */

	step = LROUND((nan - frame_len)/((float) (des_frames - 1)));

	if (step <= 0) 
	  step = 1;
      }

    /* compute (adjust, in case of -T)  number of frames to fill 
       specified segment */
    
    if (nan == 0) 
      nan = end - start + 1;

    if (nan < frame_len) 
      nframes = 1;
    else {
      nframes = 1 + (nan - frame_len) / step; 
      /* we add one more frame unless the last one ends exactly at the 
	 end point */
      if (nan > (frame_len + (nframes - 1)*step)) {
	nframes++;
	if (debug_level) 
	  Fprintf(stderr, 
	     "%s: WARNING - last frame will exceed range by %ld points\n", 
	      ProgName, frame_len + (nframes - 1)*step - nan);
      }
    }

    if (debug_level) Fprintf(stderr,
	"start = %ld, nan = %ld, frame size = %ld, step = %ld, nframes = %ld\n",
	start, nan, frame_len, step, nframes);

    /*
     * Allocate buffers
     */

    alloc_len = MAX(frame_len, tr_length);
    x = (float *) calloc(alloc_len + 1, sizeof(float));
    px = (float *) calloc(alloc_len, sizeof(float));	
    wx = (float *) calloc(alloc_len, sizeof(float));	
    y = (float *) calloc(alloc_len, sizeof(float));
    spsassert(x && y && wx, "sgram: couldn't allocate enough memory");

    if (frame_len < 1)
	ERROR_EXIT("window must be at least 1 sample wide");
    if (debug_level) Fprintf(stderr, 
	"frame length %ld, transform length %ld, no. frames %ld\n", 
	    frame_len, tr_length, nframes);

    /*
     * only read as many points as length of transform
     */

    if (tr_length < frame_len && debug_level) {
	Fprintf(stderr, 
	    "%s: WARNING - transform length %ld less than frame length %ld\n",
	    ProgName, tr_length, frame_len);
	Fprintf(stderr,"\t\t(framing determined by step size)\n");
    }	
    if (tr_length > frame_len && debug_level) {
	Fprintf(stderr, 
	    "%s: WARNING - transform length %ld greater than frame length %ld\n",
	    ProgName, tr_length, frame_len);
	Fprintf(stderr,
	    "\t %ld zeros padded for each frame\n", tr_length-frame_len);
    }

    /*
     * create, fill in, and write output file header
     */

    oh = new_header(FT_FEA);
    add_source_file(oh, iname, ih);
    oh->common.tag = YES;
    (void) strcpy(oh->common.prog, ProgName);
    (void) strcpy(oh->common.vers, Version);
    (void) strcpy(oh->common.progdate, Date);
    oh->variable.refer = iname;
    add_comment(oh, get_cmd_line(argc, argv));

    freq_format = SPFMT_SYM_EDGE;
    spec_type = SPTYP_DB;
    num_freqs = tr_length/2 + 1;

    spsassert(
	!init_feaspec_hd(oh, YES, freq_format, spec_type,
	    YES, num_freqs, SPFRM_FIXED, (float *) NULL, sf, frame_len,
	    BYTE),
	"Error filling FEA_SPEC header" )

    *add_genhd_l("start", (long *) NULL, 1, oh) = start;
    *add_genhd_l("nan", (long *) NULL, 1, oh) = nan;
    (void) add_genhd_c("sgram_method", hd_method, strlen(hd_method), oh);
    *add_genhd_l("fft_order", (long *) NULL, 1, oh) = order;
    *add_genhd_l("fft_length", (long *) NULL, 1, oh) = tr_length;
    *add_genhd_l("step", (long *) NULL, 1, oh) = step;
    *add_genhd_e("window_type", (short *) NULL, 1, window_types, oh) = win;
    *add_genhd_d("pre_emphasis", (double *) NULL, 1, oh) = preemph;
    if (Tflag) 
      *add_genhd_l("desired_frames", (long *) NULL, 1, oh) = des_frames;

    (void)update_waves_gen(ih, oh, (float) (start + frame_len/2), (float)step);

    write_header(oh, ofile);

    /*
     * set scale factor to turn into continuous spectral density
     */

    scale = 10 * log10(2.0/(tr_length * sf));

    /*
     * allocate spectral record and move to starting position
     */

    spec_rec = allo_feaspec_rec(oh, FLOAT);

    fea_skiprec(ifile, start - 1, ih);

    /*
     * main loop
     */

    position = start;
    for (i = 0; i < nframes && more; i++) {
	if (debug_level > 1)
	    Fprintf(stderr, "frame %d\n", i + 1);

	/*
	 * get sampled data and preemphasize
	 */

	if (i == 0)
	{
	    x[0] = 0.0;
	    actsize = get_sd_orecf(x + 1, (int) frame_len, (int) step, 
				   YES, ih, ifile);
	}
	else
	{
	    actsize = - 1 + get_sd_orecf(x, (int) frame_len + 1, (int) step, 
					 NO, ih, ifile);
	}
	{
	  register float *ip, *op, old, new, pre = preemph;
	  for (j = frame_len, ip = x, old = *ip++, op = px; j--; ) {
	    *op++ = (new = *ip++) - (old*pre);
	    old = new;
	  }
	}

	if ((actsize < frame_len) && debug_level) 
	  Fprintf(stderr, 
	       "%s: WARNING - got fewer points than frame_len in frame %d\n",
	       ProgName, i + 1);

	if (actsize == 0) break;

	if (debug_level > 1)
	    Fprintf(stderr, "got %d points\n", actsize);
	if (debug_level > 2)
	{
	    pr_farray("frame from get_sd_orec", frame_len + 1, x);
	    pr_farray("preemphasized frame", frame_len, px);
	}

	more = (actsize == frame_len);    	

	/*
	 * Window data.
	 */

	(void) gp_window(window_size, px, wx, win, (double *) NULL, &frame_len);

	if (tr_length > actsize) {
	    if (debug_level > 1) Fprintf(stderr, "padding zeros to frame\n");
	    for (j = actsize; j < tr_length; j++) 
		wx[j] = 0.0;
	}

	if (debug_level > 2) pr_farray("windowed frame, input to fft", 
	    MAX(frame_len, tr_length), wx);

	/*
	 * compute total power
	 */

	tot_power = 0.0;
	{
	  register float ssq, *fp;
	  register int cnt;
	  data_length = (tr_length < actsize ? tr_length : actsize);
	  for (cnt = data_length, ssq = 0.0, fp = wx; cnt-- ; fp++)
	    ssq += *fp * *fp;
	  tot_power = ssq / data_length;
	}
	if (debug_level > 1)
	    Fprintf(stderr, "total power = %g\n", tot_power);

	/*
	 * compute fft
	 */
	{
	  register float *wp = y, zero = 0.0;
	  register int cnt = tr_length;
	  do { *wp++ = zero; } while ( --cnt);
	}

	get_rfft(wx, y, order);

	if (debug_level > 2) {
	    Fprintf(stderr, "real, imag, and power from get_fft:\n");
	    for (j = 0; j < tr_length; j++) 
		Fprintf(stderr, "%f\t%f\t%g\n", wx[j], y[j], 
		    wx[j]*wx[j] + y[j]*y[j]);
	}

	/*
	 * fill in and write out spectral record
	 */

	zeros = 0;
	{
	  register float *op, *rp, *ip, zero = 0.0, rfpow, scl = scale, ten = 10.0;
	  register int cnt;
	  for (cnt = num_freqs, rp = wx, ip = y, op = spec_rec->re_spec_val;
	       cnt-- ; rp++, ip++)
	    {
	      rfpow =  (*ip * *ip) + (*rp * *rp);
	      if (rfpow > zero)
		*op++ = scl + ten * log10(rfpow);
	      else
		{
		  *op++ = scl - ten * log10(FLT_MAX);
		  zeros++;
		}
	    }
	}
	if (zeros && !zflag)
	    Fprintf(stderr,
		    "%s: WARNING - zero power at %d points in frame %ld.\n",
		    ProgName, zeros, i+1);

	if (debug_level > 2) {
	    Fprintf(stderr, "real, imag components in spec_rec:\n");
	    for (j = 0; j < num_freqs; j++) 
		Fprintf(stderr, "%f\t%f\t\n", 
		    spec_rec->re_spec_val[j], spec_rec->im_spec_val[j]);
	}
	*spec_rec->tot_power = tot_power;
	*spec_rec->tag = position;
	position += step;

	(void) put_feaspec_rec(spec_rec, oh, ofile);
    }

/*
 * put info in ESPS common
 */

    if (strcmp(oname, "<stdout>") != 0 && strcmp(iname, "<stdin>") != 0) {
	(void) putsym_s("filename", iname);
	(void) putsym_s("prog", ProgName);
	(void) putsym_i("start", (int) start);
	(void) putsym_i("nan", (int) nan);
    }    

    exit(0);
    /*NOTREACHED*/
}

/*
 * For debug printout of float arrays
 */

void pr_farray(text, n, arr)
    char    *text;
    long    n;
    float  *arr;
{
    int	    i;

    Fprintf(stderr, "%s: %s -- %d points:\n", ProgName, text, n);
    for (i = 0; i < n; i++)
    {
	Fprintf(stderr, "%f ", arr[i]);
	if (i%5 == 4 || i == n - 1) Fprintf(stderr,  "\n");
    }
}