hist.c

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

    fflush(fp);
    fclose(fp);
}

/*****************************************************************/
/* DOC-P:   half_cosine_hist()                                   */
/* DOC-PS:  store the values of cos(3pi/2) to cos(5PI/2) in hist */
half_cosine_hist(hist,num_bins,begin_bin,end_bin,height)
HIST **hist;
int num_bins, begin_bin, end_bin, 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 */
full_cosine_hist(hist,num_bins,begin_bin,end_bin,height)
HIST **hist;
int num_bins, begin_bin, end_bin, 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                      */
hist_area(hist,num_bins)
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    */
hist_character(hist,num_bins,mean,vari,area)
HIST **hist;
int num_bins, *area;
float *mean, *vari;
{
    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               */
max_hist(hist,num_bins)
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 */
average_hieght_hist(hist,num_bins,center,window)
HIST **hist;
int num_bins, center, 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     */
smooth_hist(from,to,num_bins,window)
HIST **from, **to;
int num_bins, 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(a,b)
int *a,*b;
{
  return (*a-*b);
}

median_filter(h,out,num_bins,size)  /* size must be ODD */
HIST **h,**out;                     /* h and out may be the same */
int num_bins,size;
{
  int bin,*out_vals,*temp,i,half_size,index,median;
  int comp(); /* comparison function for sorting */
  
  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(h,from,to,type)

HIST **h;
int from,to,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         */
build_normal_hist(hist,num_bins,mean,variance,total_area)
HIST **hist;
int num_bins;
float mean, variance, 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  */
hist_slope(hist,num_bins,center,factor)
HIST **hist;
int num_bins,center,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 */
do_hist(hist,num_bins,arr,arr_cnt)
HIST **hist;
int num_bins, arr_cnt;
float *arr;
{
    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++;
    }
}