hist.c

speech signal process tools

        fprintf(stderr,"Warning: unable to open %s for writing\n",fname);
        return;
    }

    for (i=0; i<num_bins; i++){
        fprintf(fp,"%5.3f %d %d\n",
		(hist1[i]->to-hist1[i]->from)/2.0+hist1[i]->from,
                hist1[i]->count, hist2[i]->count);
    }

    fflush(fp);
    fclose(fp);
}

/*****************************************************************/
/* DOC-P:   half_cosine_hist()                                   */
/* DOC-PS:  store the values of cos(3pi/2) to cos(5PI/2) in hist */
void half_cosine_hist(HIST **hist, int num_bins, int begin_bin, int end_bin, int height)
{
    int i;
    float factor, heightby2;

    factor = 1.0 / (float)(end_bin-begin_bin+1);
    heightby2 = (float)height / 2.0;
    for (i=begin_bin;(i<num_bins) && (i <= end_bin);i++){
        hist[i]->count = heightby2 * (1.0+cos((float)(i-begin_bin)*factor*M_PI+M_PI));
    }
}

/****************************************************************/
/* DOC-P:   full_cosine_hist()                                  */
/* DOC-PS:  store the values of cos(pi/2) to cos(5pi/2) in hist */
void full_cosine_hist(HIST **hist, int num_bins, int begin_bin, int end_bin, int height)
{
    int i;
    float factor, heightby2;

    factor = 1.0 / (float)(end_bin-begin_bin+1);
    heightby2 = (float)height / 2.0;
    for (i=begin_bin;(i<num_bins) && (i <= end_bin+(end_bin-begin_bin));i++){
        hist[i]->count = heightby2 * (1.0+cos((float)(i-begin_bin)*factor*M_PI+M_PI));
    }
}

/*************************************************************/
/* DOC-P:   hist_area()                                      */
/* DOC-PS:  return the area of the HIST                      */
int hist_area(HIST **hist, int num_bins)
{
    int i, sum=0;

    for (i=0; i<num_bins; i++)
        sum+=hist[i]->count;
    return(sum);
}

/*************************************************************/
/* DOC-P:   hist_character()                                 */
/* DOC-PS:  compute the mean, variance and area of a HIST    */
void hist_character(HIST **hist, int num_bins, float *mean, float *vari, int *area)
{
    int i;
    float sum=0.0, ave;

    *area = hist_area(hist,num_bins);
    for (i=0; i<num_bins; i++)
        sum+=(hist[i]->from+(hist[i]->to - hist[i]->from)/2.0) * hist[i]->count;
    *mean = sum/(float)*area;
    sum=0.0;
    for (i=0; i<num_bins; i++){
        if (hist[i]->count > 0){
            ave = hist[i]->from + (hist[i]->to - hist[i]->from)/2.0 ;
            sum+= (ave - *mean) * (ave - *mean) * hist[i]->count;
            printf(" var %f ave %f mean %f\n",sum,ave,*mean);
        }
    }
    *vari = sum / (float)*area;
}

/*************************************************************/
/* DOC-P:   max_hist()                                       */
/* DOC-PS:  return the maximum value in a HIST               */
int max_hist(HIST **hist, int num_bins)
{
    int i, max=0;

    for (i=0; i<num_bins; i++)
        if (max<hist[i]->count)
            max=hist[i]->count;
    return(max);
}

/*********************************************************************/
/* DOC-P:   average_hieght_hist()                                    */
/* DOC-PS:  return tge average hieght in a the center of this window */
int average_hieght_hist(HIST **hist, int num_bins, int center, int window)
{
    int i, sum=0;

    for (i=center-window; i<=center+window; i++)
        if ((i>=0) && (i<num_bins))
	  sum+=hist[i]->count;
    return(sum/(window*2+1));
}

/*************************************************************/
/* DOC-P:   smooth_hist()                                    */
/* DOC-PS:  compute the average hieghts of the histogram     */
void smooth_hist(HIST **from, HIST **to, int num_bins, int window)
{
    int i, value=0,window2=window*2;

    for (value=0,i=(-window); i<(num_bins+window); i++){
        if (i-window>=0)
           value -= from[i-window]->count;
        if (i+window<num_bins)
           value += from[i+window]->count;
        if ((i>=0) && (i<num_bins))
           to[i]->count = value/window2;
    }
}

int comp(const void *a, const void *b)
{
    int *ai=(int *)a, *bi=(int *)b;
    return (*ai-*bi);
}

void median_filter(HIST **h, HIST **out, int num_bins, int size)  /* size must be ODD */
                                    /* h and out may be the same */
                  
{
  int bin,*out_vals,*temp,i,half_size,index,median;
  
  half_size=size/2;

  out_vals = (int *) malloc (num_bins * sizeof (int));
  temp = (int *) malloc (size * sizeof (int));

  for (bin=0; bin <= num_bins; bin++) { /* loop through all the bins */

    /* update analysis window */
    for (i=0; i<size; i++) {
      index=bin-half_size+i;
      temp[i]=((index >= 0) && (index < num_bins)) ? h[index]->count : 0;
    }
      
    /* find median of window */
    qsort(temp,size,sizeof(int),comp);
    median = temp[half_size];

    out_vals[bin]=median;
  }

  /* write filtered values back into output histogram */
  for (i=0; i<num_bins; i++)
    out[i]->count = out_vals[i];

  free(out_vals);
  free(temp);
}

    
int locate_extremum(HIST **h, int from, int to, int type)
{
  int i,j,extremum,swing_loc,k,next_swing_loc,diff1,diff2,pre_swing,post_swing;

  for (i=from+PEAK_WIDTH; i<to-PEAK_WIDTH; i++) {
    if (h[i]->count==0) continue; /* not interested in extrema at 0 */
    extremum=1; /* assume it's an extremum to begin with */
    pre_swing=post_swing=0;
    swing_loc=i-PEAK_WIDTH;
    for (j=i-PEAK_WIDTH; j<i; j++) /* check the preceding samples */
      if (type==PEAK) {
	if (h[j]->count > h[j+1]->count) {
	  extremum=0;
	  break;
	}
	if (h[j]->count != h[j+1]->count) {
	  pre_swing=1;
	}
      }	else { /* type == TROUGH */
	if (h[j]->count < h[j+1]->count) {
	  extremum=0;
	  break;
	}
	if (h[j]->count != h[j+1]->count) {
	  pre_swing=1;
	}
      }

    if (!extremum) continue;

    for (j=i; j<i+PEAK_WIDTH; j++) /* check the subsequent samples */
      if (type==PEAK) {
	if (h[j]->count < h[j+1]->count) {
	  extremum=0;
	  break;
	}
	if (h[j]->count != h[j+1]->count) {
	  post_swing=1;
	}
      }	else { /* type == TROUGH */
	if (h[j]->count > h[j+1]->count) {
	  extremum=0;
	  break;
	}
	if (h[j]->count != h[j+1]->count) {
	  post_swing=1;
	}
      }

    /* check to make sure it isn't a step function */
    /* this kind of check is necessary if the peak is wider than the window */
    if (((pre_swing+post_swing)<=1)&&(extremum)) {
      for (k=i; k>from; k--)
	if ((diff1 = (h[k-1]->count - h[k]->count)) != 0)
	  break;
      swing_loc=k;
      for (k=i; (k<to-1); k++)  /* find next swing */
	if ((diff2 = (h[k]->count - h[k+1]->count)) != 0)
	  break;
      next_swing_loc=k;
      if ((type==PEAK)&&((diff1>0)||(diff2<0))) continue;   /* no dice */
      if ((type==TROUGH)&&((diff1<0)||(diff2>0))) continue; /* ditto   */

      /* otherwise, the peak is at the mid-point of this plateau */
      return (int) (swing_loc+next_swing_loc)/2;
    }

    if (extremum) return i;
  }

  return to;
}

/*************************************************************/
/* DOC-P:   build_normal_hist()                              */
/* DOC-PS:  store the the HIST a normal distribution         */
void build_normal_hist(HIST **hist, int num_bins, float mean, float variance, float total_area)
{
    int i;
    float x, value, con, e;

    printf("Building a histogram for a normal distribution\n");
    printf("       mean %f,  variance %f, area %f\n",mean,variance,total_area);

    con = (1.0 / sqrt(2.0*M_PI*variance));
    /* first generate a N(mean,variance) dist scaled to 10000*/
    for (i=0; i<num_bins; i++){
        x = hist[i]->from + (hist[i]->to - hist[i]->from)/2.0;
        e = exp(-( (x-mean)*(x-mean) / (2*variance)));
        value = 10000.0 * con * e;
        hist[i]->count = value;
    }
    /* now normalize it to the total_area */
    { float scale; int area;
      area = hist_area(hist,num_bins);
      scale = (float)total_area/(float)area;
      for (i=0; i<num_bins; i++)
          hist[i]->count *= scale;
    }
}

/*************************************************************/
/* DOC-P:   hist_slope()       */
/* DOC-PS:  return the slope of histogram the center  */
int hist_slope(HIST **hist, int num_bins, int center, int factor)
{
    int ind, cnt;

    for (ind=0,cnt=0; ind < factor; ind++)
        if (center-ind < 0) 
            cnt -= hist[center+ind]->count;
        else if (ind+center>=num_bins)
            cnt += hist[center-ind]->count;
        else
            cnt += hist[center-ind]->count - hist[center+ind]->count;
    return((int)-(((float)cnt/(float)factor)*1000.0));
}

/*************************************************************/
/* DOC-P:   do_hist()       */
/* DOC-PS:  using an input array, compute a historgram */
void do_hist(HIST **hist, int num_bins, float *arr, int arr_cnt)
{
    int i, index, out_of_range=0;
    float dist, from;

    from = hist[0]->from;
    dist = hist[num_bins-1]->to - hist[0]->from;
    for (i=0; i<arr_cnt; i++){
        index = (int)((float)num_bins * (*(arr+i) - from) / dist);
        if ((index>=0) && (index<num_bins)) hist[index]->count++;
        else
            out_of_range++;
    }
}