reas_spectro.c

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

  if (TFR_not_reassigned.time_instants[0] <0)
    {
      printf ("reas_spectro.c : The field TFR_not_reassigned.time_instants is not correctly set\n");
      exit(0);
    }
 /*checks the compatibility between TRF_reassigned ans TFR_not_reassigned*/

   if ((TFR_not_reassigned.N_time != TFR_reassigned.N_time ) ||
       (TFR_not_reassigned.N_freq != TFR_reassigned.N_freq ) )
    {
      printf ("reas_spectro.c : incompatible fields in TFR_reassigned and TFR_not_reassigned\n");
    }

   flag=FALSE;
   for (time=0 ; time <TFR_reassigned.N_time ; time++)
     {
       if (TFR_not_reassigned.time_instants[time]
                 != 
            TFR_reassigned.time_instants[time])
	 {
	   flag=TRUE;
	 }
     }

   if (flag == TRUE)
     {
       printf ("reas_spectro.c : incompatible fields in TFR_reassigned and TFR_not_reassigned\n");
       exit(0);
     }

  step_time = TFR_reassigned.time_instants[1] - TFR_reassigned.time_instants[0];

   flag=FALSE;
   for (time=0 ; time <(TFR_reassigned.N_time-1) ; time++)
     {
       if ((TFR_reassigned.time_instants[time+1]-
            TFR_reassigned.time_instants[time])
           != step_time)
	 {
	   flag=TRUE;
	 }
     }

  if (flag == TRUE)
     {
       printf ("reas_spectro.c : time instants not regularly spaced\n");
       exit(0);
     }

 /*                   checks that the window length is odd           */

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

 /*                     tests the weighing fiels                     */



 /*--------------------------------------------------------------------*/
 /*           creation of the vector of frequency bins  (output)       */
 /*--------------------------------------------------------------------*/
  for (time = 0; time < TFR_reassigned.N_freq; time++)
    {
      TFR_reassigned.freq_bins[time] = (double) time / TFR_reassigned.N_freq;
    }
  for (time = 0; time < TFR_not_reassigned.N_freq; time++)
    {
      TFR_not_reassigned.freq_bins[time] = (double) time / TFR_not_reassigned.N_freq;
    }

  /*-----------------------------------------------------------------*/
  /*                        MEMORY ALLOCATION                        */
  /*-----------------------------------------------------------------*/
  /* allocation of memory for the norm vectors */
  normh = (double *) ALLOC (TFR_reassigned.N_time, sizeof (double));
  normth = (double *) ALLOC (TFR_reassigned.N_time, sizeof (double));
  normdh = (double *) ALLOC (TFR_reassigned.N_time, sizeof (double));

  /* allocation of memory for the windows */
  twindow = (double *) ALLOC (Window_Length, sizeof (double));
  dwindow = (double *) ALLOC (Window_Length, sizeof (double));



  /*-----------------------------------------------------------------*/
  /*                   INITIALIZATION of the OTHER FIELDS            */
  /*-----------------------------------------------------------------*/
  stft_h.is_complex = TRUE;
  stft_th.is_complex = TRUE;
  stft_dh.is_complex = TRUE;
  stft_h.N_time = TFR_reassigned.N_time;
  stft_h.N_freq = TFR_reassigned.N_freq;
  stft_th.N_time = TFR_reassigned.N_time;
  stft_th.N_freq = TFR_reassigned.N_freq;
  stft_dh.N_time = TFR_reassigned.N_time;
  stft_dh.N_freq = TFR_reassigned.N_freq;

  /* allocation of memory for the STFTs matrices */

  /* in order to save memory, 'stft_h.real_part' is stored in the same table 
     as 'TFR_not_reassigned.real_part' */

  mem_alloc_TFR (&stft_h, NULL, TFR_reassigned.time_instants,
		 TFR_not_reassigned.real_part, NULL);
  /* in order to save memory, 'stft_th.real_part is stored in the same table 
     as 'field_time' */
  mem_alloc_TFR (&stft_th, NULL, TFR_reassigned.time_instants,
		 field_time, NULL);
  /* in order to save memory, 'stft_dh.real_part'is stored in the same table 
     as 'field_freq' */
  mem_alloc_TFR (&stft_dh, NULL, TFR_reassigned.time_instants,
		 field_freq, NULL);

  /*-----------------------------------------------------------------*/
  /*                       COMPUTATION                               */
  /*-----------------------------------------------------------------*/
  step_time = TFR_reassigned.time_instants[1] - TFR_reassigned.time_instants[0];

  /* computation of the STFT of the signal with 'window' */

  stft (Signal, Window, Window_Length, stft_h, normh);

  /* dwindow is the derivative of window */
  /* (here twindow is used as an compulsary but unusefull output) */
  gradient (Window, Window_Length, 1, 1, 1, dwindow, twindow);
  step = (Window[0] + Window[Window_Length - 1]) / 2.0;
  dwindow[0] = dwindow[0] + step;
  dwindow[Window_Length - 1] = dwindow[Window_Length - 1] - step;

  /* twindow is the window multiplied by a vector of time */
  /* -(Window_Length/2-1) .... 0 ... Window_Length/2) */
  for (time = 0; time < Window_Length; time++)
    {
      twindow[time] = Window[time] * (time - (Window_Length - 1.0) / 2.0);
    }


  /* computation of the STFT of the signal with 'twindow' */
  stft (Signal, twindow, Window_Length, stft_th, normth);

  /* computation of the STFT of the signal with 'dwindow' */
  stft (Signal, dwindow, Window_Length, stft_dh, normdh);

  /* computation of the reassignement fields */
  /* in order to save memory, the first field is stored in stft_th_real
     and the other in stft_dh_real */

  for (time = 0; time < TFR_reassigned.N_time; time++)
    {
      for (freq = 0; freq < TFR_reassigned.N_freq; freq++)
	{
	  index = idx (freq, time, TFR_reassigned.N_freq);

	  /* module = spectrogram (not reassigned) at the place 'index' */
	  module = sqr (stft_h.real_part[index]) + sqr (stft_h.imag_part[index]);
	  if  (module > EPS)
	    {
	      /* the first field equals REAL(STFT_TH / STFT_H) */
	      field_time[index] = (stft_th.real_part[index] *
				   stft_h.real_part[index]
				   + stft_th.imag_part[index] *
				   stft_h.imag_part[index])
		/ module;

	      /* the second field equals - IMAG(STFT_DH / STFT_H) */
	      field_freq[index] = -(stft_dh.imag_part[index] *
				    stft_h.real_part[index]
				    - stft_dh.real_part[index] *
				    stft_h.imag_part[index])
		/ module;


	      /* normalization of the fields of reassignement */

	      field_time[index] = field_time[index]
		* ( normh[time]/ normth[time]) / step_time;

	      field_freq[index] = field_freq[index]
		* (normh[time] / normdh[time])
		* (0.5*TFR_reassigned.N_freq / pi);

	    }
	  else
	    /* when the spectrogram is nearly null : no reassignement */
	    {
	      field_time[index] = 0.0;
	      field_freq[index] = 0.0;
	    }
	
	  /* the original (not reassigned) spectrogram is stored in stft_th_imag */
	  TFR_not_reassigned.real_part[index] = module;
 
	}
    }
  
  /* Reassignement of the spectrogram according to the field of vectors */
   reassign (TFR_not_reassigned, field_time, field_freq, TFR_reassigned); 
  

  /*-----------------------------------------------------------------*/
  /*                         FREE MEMORY                             */
  /*-----------------------------------------------------------------*/
 
   FREE (dwindow); 
   FREE (twindow); 
   FREE (normdh); 
   FREE (normth);
   FREE (normh);
   FREE (stft_dh.imag_part);
   FREE (stft_th.imag_part);
   FREE (stft_h.imag_part);
   FREE (stft_dh.freq_bins);
   FREE (stft_th.freq_bins);
   FREE (stft_h.freq_bins);
 
}