cw.c

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

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

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

	     normK=0;
	     for(mu=-mumin;mu<=mumax;mu++)
	       {
	        normK = normK +
		     CWKernel[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])
                               * CWKernel[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])
                              * CWKernel[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])
                              * CWKernel[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])
		                        * CWKernel[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])
                             * CWKernel[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])
		                       * CWKernel[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;
	     for(mu=-mumin;mu<=mumax;mu++)
        {
	      normK = normK + CWKernel[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])
		                       * CWKernel[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])
       	                    * CWKernel[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])
		                       * CWKernel[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])
		                       * CWKernel[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])
                             * CWKernel[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])
		                       * CWKernel[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 (CWKernel);
}