grd.c

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

      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 (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 +
		GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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 +
		GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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])
		                    *  GRDKernel[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 and the lacf is reinitialized */
      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 (GRDKernel);
}