bj.c

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

      time = ((int) tfr.time_instants[column]) - 1;

      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);

      if (Signal.is_complex == TRUE)
         {
	        /* the signal is complex-valued */
          lacf_real[0] = Signal.real_part[time]*Signal.real_part[time]
                       + Signal.imag_part[time]*Signal.imag_part[time];
	       lacf_imag[0]=0.0;
	      }
      else
	      {
	       /* the signal is real-valued */
	       lacf_real[0] = Signal.real_part[time]*Signal.real_part[time];
          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);
	  
	      normT=0;
	      for(row=-mumin ; row <= mumax; row++)
	        {
            normT = normT + WindowT[half_WindowT_Length + row];
	        }

         for(mu=-mumin;mu<=mumax;mu++)
	        {
	         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])
   		                       *  WindowT[half_WindowT_Length+mu]/normT;
		 
		        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])
	 	                          *  WindowT[half_WindowT_Length+mu]/normT;
		 
		        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])
                                *  WindowT[half_WindowT_Length+mu]/normT;
		 
		        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])
                                *  WindowT[half_WindowT_Length+mu]/normT;
	          }
            else
             {
		        R1_real = R1_real + (Signal.real_part[time+tau-mu]
                                *  Signal.real_part[time-tau-mu])
		                          *  WindowT[half_WindowT_Length+mu]/normT;
		 
              R1_imag = 0.0;
		 
		        R2_real = R2_real + (Signal.real_part[time-tau-mu]
                                *  Signal.real_part[time+tau-mu])
                                *  WindowT[half_WindowT_Length+mu]/normT;
		 
		        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];

        }

      /* special case of 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);
	  
	     normT=0;
        for(row=-mumin ; row <= mumax; row++)
	       {
              normT = normT + WindowT[half_WindowT_Length + row];
	       }

        for(mu=-mumin;mu<=mumax;mu++)
	       {
	        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])
   		                        *  WindowT[half_WindowT_Length+mu]/normT;
		 
		         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])
	 	                           *  WindowT[half_WindowT_Length+mu]/normT;
		 
		         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])
                                 *  WindowT[half_WindowT_Length+mu]/normT;
		 
		         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])
                                 *  WindowT[half_WindowT_Length+mu]/normT;
	          }
	     
          else
	          {
		         R1_real = R1_real + (Signal.real_part[time+tau-mu]
                                 *  Signal.real_part[time-tau-mu])
		                           *  WindowT[half_WindowT_Length+mu]/normT;
		 
		         R1_imag = 0.0;
		 
		         R2_real = R2_real + (Signal.real_part[time-tau-mu]
                                 *  Signal.real_part[time+tau-mu])
		                           *  WindowT[half_WindowT_Length+mu]/normT;
		 
		        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);

}