kernel.c

Time-Frequency Toolbox,其中包含很常用的MATLAB程序

/* EXISTS AN INTERFACE PROGRAM TO MATLAB : KERMEX.C		      *
 *====================================================================*
 * Name of the function : kernel (void)				      *
 * Author               : Manuel DAVY - IRCYN			      *
 * Date of creation     : 20 - 01 - 1999			      *
 *--------------------------------------------------------------------*
 * THE ALGORITHM          					      *
 * This function generates various kernel given a shape and the params*
 * The kernel is computed in the AMBIGUITY PLANE (doppler-delay). The *
 * kernel is a matrix of given size. The kernels are :		      *
 * - MTEK  : multiform tiltable exponential kernel 		      *
 * - RGK   : Radially gaussian kernel				      *
 * - GMCWK : Generalized marginals Choi-Williams kernel		      *
 * - WV    : Wigner-Ville kernel				      *
 * See the references given in the comments below.		      *
 *								      *
 * Dans le cadre de l'analyse temps-frequence, cette fonction genere  *
 * des noyaux de Representations temps-frequence (RTF) dans le plan   *
 * des ambiguites (retard - doppler). Le noyau est mis sous la forme  *
 * d'une matrice de dimensions passees a la fonction. Chaque noyau est*
 * defini par un type et des parametres.                              *
 *								      *
 *====================================================================*
 * INPUT VARIABLES     						      *
 * Name           |                   role                  	      *
 * N_Doppler      | Number of doppler bins i.e. number of rows in the *
 *                | kernel matrix                            	      *
 * N_Delay        | Number of delays bins i.e. number of columns in   *
 *                | the kernel matrix                        	      *
 * type           | Kernel type (see the definitions in 'tftb.h')     *
 * parameters     | Parameter vector for the kernel      )    	      *
 * nb_param       | Number of parameters                     	      *
 * ker.N_doppler  | Number of lines in the kernel matrix              *
 * ker.N_delay    | Number of columns in the kernel matrix            *
 *--------------------------------------------------------------------*
 * OUTPUT VARIABLES                                                   *
 * Name           |                   role                  	      *
 * ker            | Matrix containing the computed kernel    	      *
 *--------------------------------------------------------------------*
 * INTERNAL VARIABLES						      *
 * Name           |                   role                  	      *
 * line,col       | Line and columns index in the matrix              * 
 * doppler        | Value of the doppler (resp. delay) in the current *
 * delay          |  line (resp. column)  			      *
 * inter          | Intermediary value for a computation	      *
 *====================================================================*/



void
kernel (int type, double *parameters, int nb_param, type_AF ker)
{
  int            line, col;
  double         doppler, delay;
  double         inter;


  /* some tests */

  if ((ker.N_doppler <1) || (ker.N_delay < 1))
    {
      printf("kernel.c : invalid number of lines / columns in the kernel matrix \n");
      exit(0);
    }

  switch (type)
    {
      /***************************************************************
       *              Multiform Tiltable Exponentiel Kernel          *
       ***************************************************************
       * The parameters vector is made of                            *
       * alpha, beta, gamma, r, tau0, nu0, lambda]                   *
       *   see the reference :                                       * 
       *	H. Costa and G.F. Boudreaux-Bartels,                 *
       *	Design of Time-Frequency Representations Using a     *
       *		  Multiform, Tiltable Exponential Kernel     *
       *	IEEE Trans. on Signal Processing                     *
       *	October 1995, vol. 43, no 10, pp 2283-2301           *
       ***************************************************************/
      /* LOCAL VARIABLES                                             *
       * Name           |                   role                     *
       * A              | intermediary in the computation of the     *
       *                | MTEK kernel                                *
       ***************************************************************/
   case MTEK:
      {
	/* local variables for the MTEK */

	double         A;

	/* Kernel Construction */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    doppler = (line - ker.N_doppler / 2.0 + 1.0) / ker.N_doppler;
	    for (col = 0; col < ker.N_delay; col++)
	      {
		delay = col - ker.N_delay / 2.0 + 1.0;
		A = (doppler * delay) / (TAU_0 * NU_0);

		/* case of the symmetrical kernel */
		if ((BETA == 2) && (GAMMA == 0.5))
		  {
		    A = ABS (A);
		  }

		/* case where the marginals do not have to be verified */
		if (ALPHA == 0)
		  {
		    inter = sqr (delay / TAU_0) + sqr (doppler / NU_0)
		      + 2.0
		      * R * A;
		  }
		/* case where the marginals have to be verified */
		else
		  {
		    inter = sqr (delay / TAU_0) * powof (sqr (doppler /
							      NU_0), ALPHA)
		      + sqr (doppler / NU_0) * powof (sqr (delay /
							   TAU_0), ALPHA)
		      + 2.0 * R * A;
		  }
		/* test to avoid the computation of log(0) */
		if (inter == 0)
		  {
		    ker.real_part[idx (line, col, ker.N_doppler)] = 1.0;
		  }
		else
		  {
		    ker.real_part[idx (line, col, ker.N_doppler)] =
		      exp (-pi * powof (sqr (inter), LAMBDA));
		  }
		inter = 0;
		A = 0;
	      }
	  }
      }
      break;
      /***************************************************************
       *                 Radially Gaussian kernel                    *
       ***************************************************************
       * The parameters vector is made according to the rule :       *
       * if the order of the kernel is p, the vector is              *
       * [ c ,a1 , ... , ap, b1,... ,bp] where c is the constant     *
       * ai are the cosine coefficients and bi the sine              *
       * coefficients in the Fourier series decomposition of the     *
       * contour. See the reference :                                *
       *     M. Davy and C. Doncarli,                                *
       *    Optimal Kernels of Time-Frequency Representations for    *
       *    Signal Classification,                                   *
       *	  TFTS 1998, pp 581-584.                             *
       ***************************************************************/
      /* LOCAL VARIABLES                                             *
       * Name           |                   role                     *
       * order          | Fixes the maximum p in the vector of params*
       * p              | Current parameter p in the vector of params*
       * a,b            | Vectors containing the parameters ap, bp   *
       *                | from the vector of params                  *
       * c              | The minimum value of the contour function  *
       * inter, mini    | Intermediary values in the computations    *
       * phi            | angle parameter in the polar coordinates   *
       * rho2           | square radius parameter in polar coord.    *         
       ***************************************************************/


    case RGK:
      {
	/* local variables for the RGK */
	int            order, p;
	double         *a, *b;
	double         c, inter, mini;
	double         phi, rho2;


	/* some error cases to avoid ... */
	if (ISODD(nb_param) == 0)
	  {
	    printf("kernel.c : the number of RGK parameters must be ODD\n");
	    exit(0);
	  } 


	order = (nb_param - 1) / 2;


	/* memory allocation for a and b */
	a = (double *) ALLOC ( order , sizeof(double) );
	b = (double *) ALLOC ( order , sizeof(double) );

	/* variables recovery */
	c = parameters[0];
	for (p = 0; p < order; p++)
	  {
	    a[p] = parameters[p + 1];
	    b[p] = parameters[order + p + 1];
	  }
	/*-----------------------------------------------*/
	/*             Kernel Construction               */
	/*-----------------------------------------------*/

	/* minimum value of the contour function */
	mini = 0;	


	/* construction of the matrix of the contour function */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    /* current doppler value */
	    doppler = (line - ker.N_doppler / 2.0 + 1.0) / ker.N_doppler;
	    /* normalization of the doppler to have angles in radians */
	    doppler = doppler * sqrt (ker.N_delay);

	    for (col = 0; col < ker.N_delay; col++)
	      {
		/* currrent delay value */
		delay = col - ker.N_delay / 2.0 + 1.0;
		/* normalization of the delay to have angles in radians */
		delay = delay / sqrt (ker.N_delay);

		/* computation of the angles in the ambiguity plane */
		if (((delay > 0) && (doppler > 0))
		    || ((delay < 0) && (doppler < 0)))
		  {
		    phi = atan (doppler / delay);
		  }
		else
		  {
		    phi = atan (doppler / delay) + pi;
		  }

		inter = 0;
		for (p = 0; p < order; p++)
		  {
		    inter = inter + a[p] * cos (2.0 * (p + 1) * phi) +
		      b[p]
		      * sin (2.0 * (p + 1) * phi);
		  }
		/* look for the minimum */
		if (inter < mini)
		  {