distance.c

来自「Time-Frequency Toolbox,其中包含很常用的MATLAB程序」· C语言 代码 · 共 537 行 · 第 1/2 页

/* ***************** LQ distance ********************************** */
/*------------------------------------------------
   -                                   -  1/coef
   | / /                        coef    |
d= | | |  | TFR1(t,f)-TFR2(t,f)|    dtdf|
   | / /                                |
   -                                   -
---------------------------------------------------*/
    case LQ:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = first_TFR.real_part[index];
	      tfr2_local = second_TFR.real_part[index];
	      inter = tfr1_local - tfr2_local;
	      distan = distan + powof (ABS (inter), coef);

	    }
	}
      distan = powof (distan, 1.0 / coef);
      break;

/* ******************* Quadratic distance ************************* */
/*--------------------------------------------- 
    / /                        2    
d=  | |  | TFR1(t,f)-TFR2(t,f)| dtdf
    / /                                   
-----------------------------------------------*/
    case QUADRATIC:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = first_TFR.real_part[index];
	      tfr2_local = second_TFR.real_part[index];
	      inter = tfr1_local - tfr2_local;
	      distan = distan + inter * inter;
	    }
	}

      break;

/************************ Correlation distance ******************** */
/*---------------------------------------------  
       / /                                      
d= 1 - | | TFR1(t,f)*TFR2(t,f) dtdf
       / /   
---------------------------------------------*/
    case CORRELATION:
      distan = 0;
      first_sum = 0;
      second_sum = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index =  idx (time, freq, N_time);
	      tfr1_local = first_TFR.real_part[index];
	      tfr2_local = second_TFR.real_part[index];

	      first_sum = first_sum + tfr1_local * tfr1_local;
	      second_sum = second_sum + tfr2_local * tfr2_local;

	      inter = tfr1_local * tfr2_local;
	      distan = distan + inter;
	    }
	}
      distan = 1 - distan / (first_sum + second_sum);
      break;


/************************* Kolmogorov distance **********************/
/*--------------------------------------------- 
    / /                                    
d=  | |  | TFR1_norm(t,f)-TFR2_norm(t,f)| dtdf
    / /    
---------------------------------------------*/
    case KOLMOGOROV:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      inter = tfr1_local - tfr2_local;
	      distan = distan + ABS (inter);
	    }
	}
      break;

/*************************** Kullback distance **********************/

/*---------------------------------------------------------------
    / /                                    TFR1_norm(t,f) 
d=  | |(TFR1_norm(t,f)-TFR2_norm(t,f))*log--------------- dtdf
    / /                                    TFR2_norm(t,f)
---------------------------------------------------------------*/
    case KULLBACK:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      if ((tfr1_local != 0) && (tfr2_local != 0))
		{
		  inter = (tfr1_local - tfr2_local) * log (tfr1_local / tfr2_local);
		}
	      else
		/* the distance is not defined , then the null points do not count */
		{
		  inter = 0;
		}
	      distan = distan + ABS (inter);
	    }
	}
      break;

/*********************** Chernoff distance **************************/
/*-----------------------------------------------------------------
         | / /               coef                1/coef   |    
d= -log  | | | TFR1_norm(t,f)     *TFR2_norm(t,f)     dtdf|
         |   / /                                          |
         -                                               -
---------------------------------------------------------------*/
    case CHERNOFF:		/*Chernoff distance */
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      inter = powof (tfr1_local, coef) * powof (tfr2_local, 1.0 - coef);
	      distan = distan + inter;
	    }
	}
      distan = -log (distan);
      break;

/********************** Generalized Matusita distance ***************/
/*---------------------------------------------------------------
     -                                                       -  1/coef
    | / /               1/coef               1/coef  coef     |    
d=  | | | |TFR1_norm(t,f)      -TFR2_norm(t,f)      |     dtdf|
    | / /                                                     |
     -                                                       -
---------------------------------------------------------------*/
    case MATUSITA:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      inter = powof (tfr1_local, 1.0 / coef) - powof (tfr2_local,
							      1.0 / coef);
	      distan = distan + powof (ABS (inter), coef);
	    }
	}
      distan = powof (distan, 1.0 / coef);
      break;

/********************** Normalized Lq distance **********************/
/*---------------------------------------------------------------
   -                                             -    1/coef
   | / /                                  coef    |
d= | | |  | TFR1_norm(t,f)-TFR2_norm(t,f)|    dtdf|
   | / /                                          |
   -                                             -
---------------------------------------------------------------*/
    case NLQ:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      inter = tfr1_local - tfr2_local;
	      distan = distan + powof (ABS (inter), coef);
	    }
	}
      distan = powof (distan, 1.0 / coef);
      break;

/**************** q Normalized Log Spectral deviation ***************/
/*---------------------------------------------------------------
   -                                                       -  1/coef
   | / /                                            coef    |
d= | | |  | log(TFR1_norm(t,f))-log(TFR2_norm(t,f))|    dtdf|
   | / /                                                    |
   -                                                       -
---------------------------------------------------------------*/
    case LSD:
      distan = 0;
      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      tfr1_local = ABS (first_TFR.real_part[index]) / first_sum;
	      tfr2_local = ABS (second_TFR.real_part[index]) / second_sum;
	      inter = log (tfr1_local) - log (tfr2_local);
	      distan = distan + powof (ABS (inter), coef);
	    }
	}
      distan = powof (distan, 1.0 / coef);
      break;


/**************************** Jensen Divergence *********************/
    case JENSEN:
      /* uses Renyi and Jensen_inter_index */
      TFR_inter.N_time = N_time;
      TFR_inter.N_freq = N_freq;
      TFR_2_norm.N_time = N_time;
      TFR_2_norm.N_freq = N_freq;
      TFR_inter.is_complex = FALSE;
      TFR_2_norm.is_complex = FALSE;

      TFR_inter.real_part = (double *) ALLOC (TFR_inter.N_time *
					       TFR_inter.N_freq, sizeof (double));
      TFR_2_norm.real_part = (double *) ALLOC (TFR_2_norm.N_time *
						TFR_2_norm.N_freq, sizeof (double));


      for (time = 0; time < N_time; time++)
	{
	  for (freq = 0; freq < N_freq; freq++)
	    {
	      index = idx (freq, time, N_freq);
	      TFR_2_norm.real_part[index] = ABS (second_TFR.real_part[index])
		/ second_sum;
	      TFR_inter.real_part[index] = (ABS (first_TFR.real_part[index])
					    / first_sum
					    + TFR_2_norm.real_part[index]) / 2;
	    }
	}
      distan = Jensen_inter_index (TFR_inter, TFR_2_norm, coef);

      FREE (TFR_inter.real_part);
      FREE (TFR_2_norm.real_part);

      break;
      
    }


/*====================================================================
      the final distance is stored in the output variable
====================================================================*/
  *dist = distan;
}