ridbn.c

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

      /* maximum value of the delay in order to take the edges
	 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;
	}

   RIDBNKernel_Length = 1;
   RIDBNKernel = (double *) ALLOC( RIDBNKernel_Length,sizeof(double) );
   if (RIDBNKernel == NULL ) mexErrMsgTxt("Memory Failure");
   RIDBNKernel[0] = 1;

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

    RIDBNKernel_Length = RIDBNKernel_Length + 2;
    RIDBNKernel = (double *) realloc(RIDBNKernel, RIDBNKernel_Length*sizeof(double));
    Binomial_Kernel(RIDBNKernel, RIDBNKernel_Length);

   for(mu = -mumin ; mu <= mumax ; mu++)
    {
        normK=normK+ RIDBNKernel[tau+mu]*WindowT[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])
                  *RIDBNKernel[tau+mu]
                  *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[half_WindowT_Length + mu]/ normK;
	       }

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


    FREE(RIDBNKernel);

     tau=floor(tfr.N_freq/2);
     RIDBNKernel = (double *) ALLOC ( (2*tau+1) , sizeof(double) );
     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;

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

         RIDBNKernel[0]=1;

         normK = 0.0;
         for( p=0; p<(2*tau) ; p++ )
            {
              RIDBNKernel[p+1] = RIDBNKernel[p]*(2*tau-(p+1)+1)/(p+1);
              normK = normK + RIDBNKernel[p+1];
            }

         normK = normK +1;
         for( p=0; p<(2*tau+1) ; p++ )
            RIDBNKernel[p] = RIDBNKernel[p]/normK;

         normK = 0.0;
         for(mu = -mumin ; mu <= mumax ; mu++)
         {
          normK=normK+ RIDBNKernel[tau+mu]*WindowT[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])
                  *RIDBNKernel[tau+mu]
                  *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[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])
                   *RIDBNKernel[tau+mu]
                   *WindowT[half_WindowT_Length + mu]/ normK;
	       }

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


    }

    FREE(RIDBNKernel);

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