ridb.c

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

	{
	  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++)
	    {
	      RIDBKernel[half_WindowT_Length + mu] =
	         sqrt(1.0 - ((double)(mu*mu))/((double)(tau*tau)))
 	       * WindowT[half_WindowT_Length + mu];

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

	    }

	  /* when the kernel is nearly null */
	  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])
                                    *    RIDBKernel[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])
                                    *    RIDBKernel[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])
                                    *    RIDBKernel[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])
                                     *   RIDBKernel[half_WindowT_Length+mu]
		                     /   normK;
		}
	      /* case of real-valued signal */
	      else
		{
		  R1_real = R1_real +    (Signal.real_part[time+tau-mu]
			  	        * Signal.real_part[time-tau-mu])
		                    *     RIDBKernel[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])
		                    *    RIDBKernel[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++)
	    {
	      RIDBKernel[half_WindowT_Length + mu] =
	         sqrt(1.0 - ((double)(mu*mu))/((double)(tau*tau)))
 	       * WindowT[half_WindowT_Length + mu];

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

	    }

	  /* when the kernel is nearly null */
	  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])
                  *    RIDBKernel[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])
                   *    RIDBKernel[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])
                   *    RIDBKernel[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])
                   *   RIDBKernel[half_WindowT_Length+mu]
		                     /   normK;
		}
	      /* case of real-valued signal */
	      else
		{
		  R1_real = R1_real +    (Signal.real_part[time+tau-mu]
			  	        * Signal.real_part[time-tau-mu])
		                    *     RIDBKernel[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])
		                    *    RIDBKernel[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(RIDBKernel);
}