ridt.c

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

	 into account */
      taumax = MIN((time+half_WindowT_Length),(Signal.length-time-1+half_WindowT_Length));
      taumax = MIN(taumax,(tfr.N_freq / 2 - 1));
      taumax = MIN(taumax, half_WindowF_Length);

      /* initialization of the first local autocorrelation function */
      if (Signal.is_complex == TRUE)
	{
	  lacf_real[0] =   Signal.real_part[time] 
	                 * Signal.real_part[time]
                       +   Signal.imag_part[time]
                         * Signal.imag_part[time] 
                       *   WindowT[half_WindowT_Length];
	  /* the imag part is always zero because the imag part
	     of any complex 'x * conjugate(x)' is zero */
	    
	  lacf_imag[0] =  0.0 ;
	}
      else /* the signal is real-valued */
	{
	  lacf_real[0] =   Signal.real_part[time] 
                         * Signal.real_part[time]
                       *   WindowT[half_WindowT_Length];

	  lacf_imag[0] =   0.0;
	}

      /* The signal is windowed around the current time */
      for (tau = 1; tau <= taumax; tau++)
	{
         R1_real=0.0;
	 R2_real=0.0;
         R1_imag=0.0;
	 R2_imag=0.0;

	 /* bound of mu in order to take into account the edges */
	 mumin=MIN(tau,half_WindowT_Length);
         mumin=MIN(mumin,(Signal.length-time-1-tau));

         mumax=MIN(tau,half_WindowT_Length);
         mumax=MIN(mumax,time-tau);


         normK = 0.0;
	 for(mu = -mumin ; mu <= mumax ; mu++)       
	   { 
	    RIDTKernel[half_WindowT_Length + mu] = 
	      (1.0 - ABS((double)mu/(double)tau))
	      *WindowT[half_WindowT_Length + mu];

	     /* computation of the current kernel norm */
	     normK = normK +  RIDTKernel[half_WindowT_Length + mu];
           }

	 /* when the kernel is nearly null, no normalization */   
         if (normK<EPS)
	   normK=1.0;
            
	 for(mu = -mumin ; mu <= mumax ; mu++)
	   {
	     /* case of complex valued signal */
	     if (Signal.is_complex == TRUE)
	       {
		  R1_real = R1_real +   (Signal.real_part[time+tau-mu]
			 	      *  Signal.real_part[time-tau-mu]
				     +   Signal.imag_part[time+tau-mu]
				      *  Signal.imag_part[time-tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
		                    /    normK;

		  R1_imag = R1_imag +   (Signal.imag_part[time+tau-mu]
				       * Signal.real_part[time-tau-mu]
				     -   Signal.real_part[time+tau-mu]
				       * Signal.imag_part[time-tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
                                    /    normK;

		  R2_real = R2_real +   (Signal.real_part[time-tau-mu]
				       * Signal.real_part[time+tau-mu]
				     +   Signal.imag_part[time-tau-mu]
				       * Signal.imag_part[time+tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
		                    /    normK;

		  R2_imag = R2_imag +   (Signal.imag_part[time-tau-mu]
				       * Signal.real_part[time+tau-mu]
                                     -   Signal.real_part[time-tau-mu]
				       * Signal.imag_part[time+tau-mu])
                                     *   RIDTKernel[half_WindowT_Length+mu]
		                     /   normK;
	       }
	     else
	       {
		 R1_real = R1_real +    (Signal.real_part[time+tau-mu]
			  	       * Signal.real_part[time-tau-mu])
		                    *    RIDTKernel[half_WindowT_Length+mu]
                                   /     normK;
		 R1_imag = 0.0;

		 R2_real = R2_real +   (Signal.real_part[time-tau-mu]
				      * Signal.real_part[time+tau-mu])
		                   *    RIDTKernel[half_WindowT_Length+mu]
		                   /    normK;
		 R2_imag = 0.0;
	       }

	   }
	  lacf_real[tau] = R1_real * WindowF[half_WindowF_Length+tau];
	  lacf_imag[tau] = R1_imag * WindowF[half_WindowF_Length+tau];
	  lacf_real[tfr.N_freq-tau] = R2_real
	                            * WindowF[half_WindowF_Length-tau];
	  lacf_imag[tfr.N_freq-tau] = R2_imag
	                            * WindowF[half_WindowF_Length-tau];
        }


     tau=floor(tfr.N_freq/2);
     if ((time<=Signal.length-tau-1)&(time>=tau)&(tau<=half_WindowF_Length))
       {
       R1_real=0.0;
	    R2_real=0.0;
       R1_imag=0.0;
	    R2_imag=0.0;

	 /* bound of mu in order to take into account the edges */
	 mumin=MIN(tau,half_WindowT_Length);
         mumin=MIN(mumin,(Signal.length-time-1-tau));

         mumax=MIN(tau,half_WindowT_Length);
         mumax=MIN(mumax,time-tau);


         normK = 0.0;
	 for(mu = -mumin ; mu <= mumax ; mu++)       
	   { 
	    RIDTKernel[half_WindowT_Length + mu] = 
	      (1.0 - ABS((double)mu/(double)tau))
	      *WindowT[half_WindowT_Length + mu];

	     /* computation of the current kernel norm */
	     normK = normK +  RIDTKernel[half_WindowT_Length + mu];
           }

	 /* when the kernel is nearly null, no normalization */   
         if (normK<EPS)
	   normK=1.0;
            
	 for(mu = -mumin ; mu <= mumax ; mu++)
	   {
	     /* case of complex valued signal */
	     if (Signal.is_complex == TRUE)
	       {
		  R1_real = R1_real +   (Signal.real_part[time+tau-mu]
			 	      *  Signal.real_part[time-tau-mu]
				     +   Signal.imag_part[time+tau-mu]
				      *  Signal.imag_part[time-tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
		                    /    normK;

		  R1_imag = R1_imag +   (Signal.imag_part[time+tau-mu]
				       * Signal.real_part[time-tau-mu]
				     -   Signal.real_part[time+tau-mu]
				       * Signal.imag_part[time-tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
                                    /    normK;

		  R2_real = R2_real +   (Signal.real_part[time-tau-mu]
				       * Signal.real_part[time+tau-mu]
				     +   Signal.imag_part[time-tau-mu]
				       * Signal.imag_part[time+tau-mu])
                                    *    RIDTKernel[half_WindowT_Length+mu]
		                    /    normK;

		  R2_imag = R2_imag +   (Signal.imag_part[time-tau-mu]
				       * Signal.real_part[time+tau-mu]
                                     -   Signal.real_part[time-tau-mu]
				       * Signal.imag_part[time+tau-mu])
                                     *   RIDTKernel[half_WindowT_Length+mu]
		                     /   normK;
	       }
	     else
	       {
		 R1_real = R1_real +    (Signal.real_part[time+tau-mu]
			  	       * Signal.real_part[time-tau-mu])
		                    *    RIDTKernel[half_WindowT_Length+mu]
                                   /     normK;
		 R1_imag = 0.0;

		 R2_real = R2_real +   (Signal.real_part[time-tau-mu]
				      * Signal.real_part[time+tau-mu])
		                   *    RIDTKernel[half_WindowT_Length+mu]
		                   /    normK;
		 R2_imag = 0.0;
	       }

	   }
	  lacf_real[tau] = 0.5*(R1_real * WindowF[half_WindowF_Length+tau]
                          +R2_real * WindowF[half_WindowF_Length-tau]);
	  lacf_imag[tau] = 0.5*(R1_imag * WindowF[half_WindowF_Length+tau]
                         +R2_imag * WindowF[half_WindowF_Length-tau]);
	  }


      /* fft of the local autocorrelation function lacf */
      fft (tfr.N_freq, Nfft, lacf_real, lacf_imag);

      /* the fft is put in the tfr matrix  */
      for (row = 0; row < tfr.N_freq; row++)
	{
	  tfr.real_part[idx (row,column,tfr.N_freq)]= lacf_real[row];
	  lacf_real[row] = 0.0;
	  lacf_imag[row] = 0.0;
        }
    }
 /*--------------------------------------------------------------------*/
 /*                free the memory used in this program                */
 /*--------------------------------------------------------------------*/
  FREE (lacf_real);
  FREE (lacf_imag);
  FREE(RIDTKernel);
}