bud.c

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

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

	      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(half_WindowT_Length, (Signal.length-time-1-tau));
	      mumax = MIN(half_WindowT_Length,time-tau);

         normK = 0.0;
         for(mu=-mumin;mu<=mumax;mu++)
           {
	         normK = normK + BUDKernel[idx(tau-1,half_WindowT_Length+mu,
			                           MIN(tfr.N_freq / 2,half_WindowF_Length))];
	        }

         for(mu=-mumin;mu<=mumax;mu++)
	        {
	        /* case of complex valued signal */

	         if (Signal.is_complex == TRUE)
             {
		         index = idx (tau - 1, half_WindowT_Length + mu,
			      MIN (tfr.N_freq / 2, half_WindowF_Length));

		         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])
		                           *  BUDKernel[index]/ 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])
		                           *  BUDKernel[index] / normK;


		         index = idx (tau - 1, half_WindowT_Length - mu,
			      MIN(tfr.N_freq / 2, half_WindowF_Length));

		         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])
		                           *  BUDKernel[index] / 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])
		                           *  BUDKernel[index] / normK;

	          }
	           /* case of real-valued signal */
	         else
	       {
		      index = idx (tau - 1, half_WindowT_Length + mu,
			                        MIN (tfr.N_freq / 2, half_WindowF_Length));


		      R1_real = R1_real + (Signal.real_part[time + tau - mu]
                              *  Signal.real_part[time - tau - mu])
		                        *  BUDKernel[index] / normK;

		      R1_imag = 0.0;


		      index = idx (tau - 1, half_WindowT_Length - mu,
			                         MIN(tfr.N_freq / 2, half_WindowF_Length));

		      R2_real = R2_real + (Signal.real_part[time - tau - mu] *
			                        Signal.real_part[time + tau - mu])
		                        *  BUDKernel[index]/ 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(half_WindowT_Length, (Signal.length-time-1-tau));
	    mumax = MIN(half_WindowT_Length,time-tau);

       normK = 0.0;
       for(mu=-mumin;mu<=mumax;mu++)
	     {
	      normK = normK +
         BUDKernel[idx(tau-1,half_WindowT_Length+mu,
                            MIN(tfr.N_freq / 2,half_WindowF_Length))];
	    }

       for(mu=-mumin;mu<=mumax;mu++)
	    {
	     /* case of complex valued signal */

	     if (Signal.is_complex == TRUE)
	       {
		      index = idx (tau - 1, half_WindowT_Length + mu,
			                   MIN (tfr.N_freq / 2, half_WindowF_Length));

		      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])
		                         *  BUDKernel[index]/ 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])
		                        *  BUDKernel[index] / normK;


		      index = idx (tau - 1, half_WindowT_Length - mu,
			                         MIN(tfr.N_freq / 2, half_WindowF_Length));

            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])
		                        *  BUDKernel[index] / 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])
		                      *  BUDKernel[index] / normK;

	       }
	     /* case of real-valued signal */
	     else
	       {
		       index = idx (tau - 1, half_WindowT_Length + mu,
			      MIN (tfr.N_freq / 2, half_WindowF_Length));


		    R1_real = R1_real + (Signal.real_part[time + tau - mu]
                            * Signal.real_part[time - tau - mu])
		                      * BUDKernel[index] / normK;

		    R1_imag = 0.0;


		    index = idx (tau - 1, half_WindowT_Length - mu,
			      MIN(tfr.N_freq / 2, half_WindowF_Length));

		    R2_real = R2_real + (Signal.real_part[time - tau - mu]
                            * Signal.real_part[time + tau - mu])
		                      * BUDKernel[index]/ 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 (BUDKernel);

}