source: trunk/SRC/ToBeReviewed/STATISTICS/c_timecorrelate.pro @ 177

;+
; @file_comments
; 
;
; @categories
; Statistics
;
; @param XD
;
;
; @param YD
;
;
; @param M
;
;
; @param NT
;
;
; @param NDIM
;
;
; @keyword ZERO2NAN
;
;
; @keyword DOUBLE
; If set to a non-zero value, computations are done in
; double precision arithmetic.
;
; @examples
;
;
; @history
;
;
; @version
; $Id$
;
;-
FUNCTION timecross_cov, Xd, Yd, M, nT, Ndim, Double = Double, ZERO2NAN = zero2nan
  ;Sample cross covariance function.

  compile_opt hidden
;
   case Ndim OF
      1:res = TOTAL(Xd[0:nT - M - 1L] * Yd[M:nT - 1L] $
                    , Double = Double)
      2:res = TOTAL(Xd[*, 0:nT - M - 1L] * Yd[*, M:nT - 1L] $
                    , Ndim, Double = Double)
      3:res = TOTAL(Xd[*, *, 0:nT - M - 1L] * Yd[*, *, M:nT - 1L] $
                    , Ndim, Double = Double)
      4:res = TOTAL(Xd[*, *, *, 0:nT - M - 1L] * Yd[*, *, *, M:nT - 1L] $
                    , Ndim, Double = Double)
   ENDCASE
   if keyword_set(zero2nan) then begin
      zero = where(res EQ 0)
      if zero[0] NE -1 then res[zero] = !values.f_nan
   ENDIF
;
   RETURN, res

END
;+
; @file_comments
; This function computes the "time cross correlation" Pxy(L) or
; the "time cross covariance" between 2 arrays (this is some
; kind of c_correlate but for multidimenstionals arrays) as a
; function of the lag (L).
;
; @categories
; Statistics
;
; @param X {in}{required}{type=array}
; An Array which last dimension is the time dimension of
; size n, float or double.
;
; @param Y {in}{required}{type=array}
; An Array which last dimension is the time dimension of
; size n, float or double.
;
; @param LAG {in}{required}{type=scalar or vector}
; A scalar or n-element vector, in the interval [-(n-2), (n-2)],
; of type integer that specifies the absolute distance(s) between
; indexed elements of X.
;
; @keyword COVARIANCE
; If set to a non-zero value, the sample cross 
; covariance is computed.
;
; @keyword DOUBLE
; If set to a non-zero value, computations are done in
; double precision arithmetic.
;
; @examples
; 
;       Define two n-element sample populations.
;         x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
;         y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95]
;
;       Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7
;         lag = [-5, 0, 1, 5, 6, 7]
;         result = c_timecorrelate(x, y, lag)
;
;       The result should be:
;         [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
;
; @history
;       - 01/03/2000 Sebastien Masson (smasson\@lodyc.jussieu.fr)
;       Based on the C_CORRELATE procedure of IDL
;       - August 2003 Sebastien Masson 
;       update according to the update made in C_CORRELATE by
;       W. Biagiotti and available in IDL 5.5
;
;       INTRODUCTION TO STATISTICAL TIME SERIES
;       Wayne A. Fuller
;       ISBN 0-471-28715-6
;
; @version
; $Id$
;
;-
FUNCTION c_timecorrelate, X, Y, Lag, Covariance = Covariance, Double = Double

;Compute the sample cross correlation or cross covariance of 
;(Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l).

   ON_ERROR, 2

   xsize = SIZE(X)
   ysize = SIZE(Y)
   nt = float(xsize[xsize[0]])
   NDim = xsize[0]

   if total(xsize[0:xsize[0]] NE ysize[0:ysize[0]]) NE 0 then $
    MESSAGE, "X and Y arrays must have the same size and the same dimensions"

;Check length.
   if nt lt 2 then $
    MESSAGE, "Time dimension of X and Y arrays must contain 2 or more elements."
   
;If the DOUBLE keyword is not set then the internal precision and
;result are identical to the type of input.
   if N_ELEMENTS(Double) eq 0 then $
    Double = (Xsize[Xsize[0]+1] eq 5 or ysize[ysize[0]+1] eq 5)

   if n_elements(lag) EQ 0 then lag = 0
   nLag = N_ELEMENTS(Lag)

;Deviations
   if double then one = 1.0d ELSE one = 1.0
   Ndim = size(X, /n_dimensions)
   Xd = TOTAL(X, Ndim, Double = Double) / nT
   Xd = X - Xd[*]#replicate(one,  nT)
   Yd = TOTAL(Y, Ndim, Double = Double) / nT
   Yd = Y - Yd[*]#replicate(one,  nT)

   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.

   case NDim of
      1:if Double eq 0 then  Cross = FLTARR(nLag) else  Cross = DBLARR(nLag)
      2:if Double eq 0 then  Cross = FLTARR(Xsize[1], nLag) else  Cross = DBLARR(Xsize[1], nLag)
      3:if Double eq 0 then  Cross = FLTARR(Xsize[1], Xsize[2], nLag) $
      else  Cross = DBLARR(Xsize[1], Xsize[2], nLag)
      4:if Double eq 0 then  Cross = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $
      else  Cross = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag)
   endcase

   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross  Crossation.
      for k = 0, nLag-1 do begin
         if Lag[k] ge 0 then BEGIN 
            case NDim of
               1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
               2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
               3: Cross[*, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
               4: Cross[*, *, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) 
             endcase
         ENDIF else BEGIN 
            case NDim of
               1: Cross[k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
               2: Cross[*, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
               3: Cross[*, *, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
               4: Cross[*, *, *, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
             endcase
         ENDELSE 
       ENDFOR
       div = sqrt(TimeCross_Cov(Xd, Xd, 0L, nT, Ndim, Double = Double, /zero2nan) * $
                  TimeCross_Cov(Yd, Yd, 0L, nT, Ndim, Double = Double, /zero2nan))
       Cross = temporary(Cross)/((temporary(div))[*]#replicate(one, nLag))
   endif else begin             ;Compute Cross Covariance.
      for k = 0, nLag-1 do begin
         if Lag[k] ge 0 then BEGIN
            case NDim of
               1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
               2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
               3: Cross[*, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
               4: Cross[*, *, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
            ENDCASE
         ENDIF else BEGIN 
            case NDim of
               1: Cross[k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
               2: Cross[*, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
               3: Cross[*, *, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
               4: Cross[*, *, *, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
            ENDCASE
         ENDELSE 
      endfor
   endelse

   if Double eq 0 then RETURN, FLOAT(Cross) else RETURN,  Cross

END