sp.c

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

/* EXISTS AN INTERFACE PROGRAM TO MATLAB : SPMEX.C                    *
 *====================================================================*
 * Name of the function : sp.c (void)                                 *
 * Authors              : Emmanuel ROY - Manuel DAVY                  *
 * Date of creation     : 10 - 02 - 1999                              *
 *--------------------------------------------------------------------*
 * THE ALGORITHM                                                      *
 *                                                                    *
 * Given a signal to analyze in time and frequency, computes the      *
 * Spectrogram (SP) :                                                 *
 *                                                                    *
 *              | /             -j 2 pi f s      |2                   *
 *    SP(t,f) = | | x(s)h*(s-t)e             ds  |                    *
 *              | /                              |                    *
 *                                                                    *
 * This function is real valued. Its computation requires a window,   *
 * its displacement positions and the number of frequency bins to be  *
 * computed with discrete variables.                                  *
 *                                                                    *
 *====================================================================*
 * INPUT VARIABLES                                                    *
 * Name              |                role                            *
 * Signal            |   The signal to analyze. No field modified     *
 *                   |                                                *
 * Window            |   Vector containing the points of the window   *
 * Window_Length     |   Number of points of the window (ODD number !)*
 *                   |                                                *
 * tfr               |  Matrix containing the resulting TFR (real)    *
 * tfr.time_instants |  positions of the smoothing window             *   
 * tfr.N_time        |  length of '.time_instants' = number of cols.  *
 *                   |  in the tfr matrix                             *
 * tfr.N_freq        |  number of frequency bins = number of rows     *
 *                   |  in the tfr matrix                             *
 * tfr.is_complex    |  must be set to FALSE (a sp is real !)         *
 *--------------------------------------------------------------------*
 * OUTPUT VARIABLES                                                   *
 * Name              |                role                            *
 * tfr.real_part     | The output tfr matrix  (real_part)             *
 * tfr.freq_bins     | Vector of normalized frequencies where the     *
 *                   | tfr is computed                                *
 * norm_vector       | Value of the normalization factor applied at   *
 *                   |  the points of computation of the sp i.e.      *
 *                   |  tfr.time_instants                             *
 *--------------------------------------------------------------------*
 * INTERNAL VARIABLES                                                 *
 * Name              |                 role                           *
 *                   |                                                *
 * Nfft              | Next power of two to tfr.N_freq                *
 *                   |                                                *
 * column, row       | variables of displacement in the matrices      *
 * tau               | local time variable (the 's' in the equation   *
 *                   | above                                          *
 * time              | Current instant of computation of the spectro  *
 * taumin            | local time variable bounds. Used to take into  *
 * taumax            | accound the beginning and the end of the       *
 *                   | signal, where the window is cut                *
 *                   |                                                *
 * half_Window_Length| half of the frequency smoothing window         *
 * normh             | current normalization value : depends on       *
 *                   | wether the window is cut or not (near the      *
 *                   | edges)                                         *
 * wind_sig_real     | Real and imaginary parts of the windowed       *
 * wind_sig_imag     | signal at the current position of the window   *
 * inter*            | several intermediary variables                 *
 *====================================================================*
 * SUBROUTINES USED HERE                                              *
 *--------------------------------------------------------------------*
 * Name   | int idx(int line, int row, int nb_row)                    *
 * Action | computes the vector index for an element in a matrix given*
 *        | the line and column indices and the number of lines       *
 * Place  | divers.c                                                  *
 *--------------------------------------------------------------------*
 * Name   | int irem( double x, double y)                             *
 * Action | computes the remainder after Euclidean division of double *
 * Place  | divers.c                                                  *
 *--------------------------------------------------------------------*
 * Name   | void fft(int n, int m, double *x, double *y)              *
 * Action | Computes the fft                                          *
 * Place  | divers.c                                                  *
 *--------------------------------------------------------------------*
 * Name   | int po2(int x)                                            *
 * Action | Computes the next power of two of x                       *
 * Place  | divers.c                                                  *
 *====================================================================*/

void
sp (type_signal Signal,
      double *Window, int Window_Length,
      type_TFR tfr, double *norm_vector)

{
  int            Nfft, column, row, time;
  int            taumin, taumax, tau;
  int            half_Window_Length;
  double        *wind_sig_real, *wind_sig_imag; /* windowed signal */
  double         inter_real, inter_imag, normh;
  double         inter;

 /*--------------------------------------------------------------------*/
 /*                      Test the input variables                      */
 /*--------------------------------------------------------------------*/


   if (tfr.is_complex == TRUE)
    {
      printf ("sp.c : The tfr matrix must be real valued\n");
      exit(0);
    }

  if (tfr.N_freq <= 0)
    {
      printf ("sp.c : The field tfr.N_freq is not correctly set\n");
      exit(0);
    }

  if (tfr.N_time <= 0)
    {
      printf ("sp.c : The field tfr.N_time is not correctly set\n");
      exit(0);
    }

 /*--------------------------------------------------------------------*/
 /*                   checks that the window length is odd             */
 /*--------------------------------------------------------------------*/

  if (ISODD(Window_Length) == 0)
    {
      printf ("sp.c : The window Length must be an ODD number\n");
      exit(0);
    }

  half_Window_Length = (Window_Length - 1) / 2;
 /*--------------------------------------------------------------------*/
 /*           creation of the vector of frequency bins  (output)       */
 /*--------------------------------------------------------------------*/
  Nfft = po2 (tfr.N_freq);

  for (row = 0; row < tfr.N_freq; row++)
    {
      tfr.freq_bins[row] = (double) row / tfr.N_freq;
    }
 /*--------------------------------------------------------------------*/
 /*                memory allocation for the windowed signal           */
 /*--------------------------------------------------------------------*/
  wind_sig_real = (double *) ALLOC (tfr.N_freq, sizeof (double));
  wind_sig_imag = (double *) ALLOC (tfr.N_freq, sizeof (double));

 /*--------------------------------------------------------------------*/
 /*      computation of the fft for the current windowed signal        */
 /*--------------------------------------------------------------------*/
  for (column = 0; column < tfr.N_time; column++)
    {
      /* initialization of the intermediary vectors */
      for (row = 0; row < tfr.N_freq; row++)
	{
	  wind_sig_real[row] = 0.0;
	  wind_sig_imag[row] = 0.0;
        }

      /* time instants of interest to compute the sp */
      time = ((int) tfr.time_instants[column]) - 1;

      /* the signal is multipied by the window between the instants
         time-taumin and time+taumax */
      /* when the window is wider than the number of desired frequencies (tfr.N_freq),
         the range is limited to tfr.N_freq */
      taumin = MIN (tfr.N_freq / 2, half_Window_Length);
      taumin = MIN (taumin, time);

      taumax = MIN ((tfr.N_freq / 2 - 1), half_Window_Length);
      taumax = MIN (taumax, (Signal.length - time - 1));

      /* Computation of a normalization factor, 
         equal to the quadratic norm of the window */

      inter = 0.0;
      for (row = -taumin; row <= taumax; row++)
	{
	  inter = inter + sqr (Window[half_Window_Length + row]);
	}
      normh = sqrt (inter);
      norm_vector[column] = 1.0 / normh;
      /* The signal is windowed around the current time */
      for (tau = -taumin; tau <= taumax; tau++)
	     {
	      row = irem( (tfr.N_freq+tau), tfr.N_freq ) ;
	      wind_sig_real[row] = Signal.real_part[time + tau]
	                         * Window[half_Window_Length + tau] 
                                 / normh;
	      if (Signal.is_complex == TRUE)
	       {
	        wind_sig_imag[row] = Signal.imag_part[time + tau]
                                   * Window[half_Window_Length + tau] 
                                   / normh;
	       }
        }
      
      /* fft of the windowed signal */
      fft (tfr.N_freq, Nfft, wind_sig_real, wind_sig_imag);


      /* the first half of the fft is put in the sp matrix  */
      for (row = 0; row < tfr.N_freq; row++)
	{
	  tfr.real_part[idx (row,column,tfr.N_freq)] = 
	             wind_sig_real[row]*wind_sig_real[row]
                   + wind_sig_imag[row]*wind_sig_imag[row];
        }
    }
 /*--------------------------------------------------------------------*/
 /*                free the memory used in this program                */
 /*--------------------------------------------------------------------*/
  FREE (wind_sig_real);
  FREE (wind_sig_imag);

}