zam.c

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

 /*--------------------------------------------------------------------*/
  for (column = 0; column < tfr.N_time; column++)
    {

      /* time instants of interest to compute the tfr */
      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);

     /* 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];
	  /* 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];

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

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

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

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

		 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];
	       }
	     else
	       {
		 R1_real = R1_real +   (Signal.real_part[time+tau-mu]
				      * Signal.real_part[time-tau-mu])
                                     *  WindowT[half_WindowT_Length+mu];

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

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

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

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

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

		 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];
	       }
	     else
	       {
		     R1_real = R1_real +   (Signal.real_part[time+tau-mu]
				      * Signal.real_part[time-tau-mu])
                                     *  WindowT[half_WindowT_Length+mu];

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

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

}