kernel.c

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

		    mini = inter;
		  }
		/* matrix of the contour function : each element in this */
		/* matrix contains the value of the contour function for */
		/* the corresponding delay and doppler values */
		ker.real_part[idx (line, col, ker.N_doppler)] = inter;

	      }
	  }

	/* construction of the RGK matrix */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    /* current normalized doppler */
	    doppler = (line - ker.N_doppler / 2.0 + 1.0) / ker.N_doppler;
	    doppler = doppler * sqrt (ker.N_delay);

	    for (col = 0; col < ker.N_delay; col++)
	      {
		/* current normalized delay */
		delay = col - ker.N_delay / 2.0 + 1.0;
		delay = delay / sqrt (ker.N_delay);
		/* Square Polar radius */
		rho2 = sqr (doppler) + sqr (delay);
		ker.real_part[idx (line, col, ker.N_doppler)] =
		  exp (-(rho2 / 2.0) / 
		       sqr (ker.real_part[idx (line, col, ker.N_doppler)]
			    - mini + c));

		/* case of the center of the ambiguity plane */
		if ((delay == 0) && (doppler == 0))
		  ker.real_part[idx (line, col, ker.N_doppler)] = 1.0;
	      }
	  }
	/* free the memory used here */
	FREE (a);
	FREE (b);
      }
      break;

      /***************************************************************
       *          Generalized Marginals Choi-Williams                *
       ***************************************************************
       * The parameters vector is made of the width of the           *
       * branches "sigma" and the angle of each branch theta_i       *
       * (0 to Pi). The number of branches is given by the number    *
       * of coefficients                                             *
       * [ sigma theta_1 theta_2  ... theta_p]                       *
       *                                                             *
       * see the reference :                                         *
       * X.-G. Xia and Y. Owechko and B. H. Soffer and R. M. Matic,  *
       * Generalized-Marginal Time-Frequency Distributions,          *
       * TFTS 1996, pp. 509-51                                       *
       ***************************************************************/
      /* LOCAL VARIABLES                                             *
       * Name           |                   role                     *
       * N_Branch       | Number of kernel branches                  *
       * branch         | Current branch number                      *
       * angle          | vector of branches angles                  *
       * sigma          | value of the parameter sigma (kernel width)*
       * inter          | intermediary variable                      *
       ***************************************************************/


    case GMCWK:
      {
	/* local variables for the RGK */
	int            N_branch, branch;
	double         *angle, sigma, inter;



	/* some error cases to avoid ... */
	if (nb_param < 2)
	  {
	    printf("kernel.c : at least 2 GMCWK parameters required\n");
	    exit(0);
	  } 


	/* variables recovery */
	sigma = parameters[0];
	N_branch = nb_param - 1;



	/* recovery of the angles */
	angle = &(parameters[1]);

	/* Kernel Construction  */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    doppler = (line - ker.N_doppler / 2.0 + 1.0) 
	      / ker.N_doppler;
	    doppler = doppler * sqrt (ker.N_delay);

	    for (col = 0; col < ker.N_delay; col++)
	      {
		delay = col - ker.N_delay / 2.0 + 1.0;
		delay = delay / sqrt (ker.N_delay);
		inter = 1;
		for (branch = 0; branch < N_branch; branch++)
		  {
		    inter = inter 
		      * sqr (doppler * cos (angle[branch])
			     + delay * sin (angle[branch]));
		  }
		ker.real_part[idx (line, col, ker.N_doppler)] = 
		  exp (-inter / sigma);
	      }
	  }
      }
      break;
      /***************************************************************
       *                Wigner-Ville kernel                          *
       ***************************************************************
       * No parameter is required. The matrix is equal to            *
       * one in each point                                           *
       ***************************************************************/


    case WIGNER:
      {



	/* Kernel Construction  */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    for (col = 0; col < ker.N_delay; col++)
	      {
		ker.real_part[idx (line, col, ker.N_doppler)] = 1.0;
	      }
	  }
      }
      break;
      /***************************************************************
       *                   Spectrogram  kernel                       *
       ***************************************************************
       *  The paremters are the window points in the time            *
       * domain. The kernel is the ambiguity function of the window  *
       ***************************************************************/
      /* LOCAL VARIABLES                                             *
       * Name           |                   role                     *
       * window         | structure containing the window of the     *
       *                | spectrogram                                *
       * AF_h           | Ambiguity function of th window            *
       * min_delay |    | in order to pad with zeros, defines the    *
       * max_delay |    | indices where the zero padding begins/ends *
       * index          | indice to access the elements in matrices  *
       *                | (stored as vectors)                        *
       ***************************************************************/



    case SPECTRO:
      {
	/* local variable */
	type_signal    window;
	type_AF        AF_h;
	int            min_delay, max_delay, index;

	/* some error cases to avoid ... */
	if (ISODD(nb_param) == 1)
	  {
	    printf("kernel.c : the window length must be EVEN for SPECTRO\n");
	    exit(0);
	  } 



	/* initialization of the local variables : window */
	window.length = ker.N_delay;
	window.is_complex = TRUE;
	/* creation of memory */
	window.real_part = (double *) ALLOC (window.length, sizeof (double));
	window.imag_part = (double *) ALLOC (window.length, sizeof (double));


	min_delay = (int) (window.length / 2 - nb_param / 2);
	max_delay = (int) (window.length / 2 + nb_param / 2 - 1);

	/* initialization of the imaginary part for the window 
	   and zero padding out of [min_delay;max_delay] */
	for (col = 0; col < window.length; col++)
	  {
	    if ((col >= min_delay) && (col <= max_delay))
	      {
		window.real_part[col] = parameters[col - min_delay];
	      }
	    else
	      {
		window.real_part[col] = 0;
	      }
	    window.imag_part[col] = 0;
	  }


	/* initialization of the local variables : AF_h */
	AF_h.N_doppler = ker.N_doppler;
	AF_h.N_delay = ker.N_delay;
	AF_h.is_complex = TRUE;

	mem_alloc_AF (&AF_h, NULL, NULL, NULL, NULL);
 	for (index = 0; index < AF_h.N_delay; index++)
	  {
	    AF_h.delay_bins[index] = -AF_h.N_delay / 2.0 + 1.0 + index;
	  }




	/* computation of the AF of the window */
	af (window, AF_h);

	/* Warning : if ker.N_delay != ker.N_doppler,
	   one must add or suppress lines in this AF */
	for (line = 0; line < ker.N_doppler; line++)
	  {
	    for (col = 0; col < ker.N_delay; col++)
	      {
		index = idx (line, col, ker.N_delay);
		ker.real_part[idx (line, col, ker.N_doppler)] =
		  sqrt (sqr (AF_h.real_part[index])
			+ sqr (AF_h.imag_part[index]));
	      }
	  }
	/* free memory !! */

	mem_free_AF (&AF_h);
	FREE (window.real_part);
	FREE (window.imag_part);

      }
      break;
    }
}