source: trunk/ToBeReviewed/STATISTICS/a_correlate2d.pro @ 21

Last change on this file since 21 was 21, checked in by pinsard, 18 years ago

upgrade of STATISTICS according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 3.5 KB
Line 
1;+
2; NAME:
3;       A_CORRELATE2d
4;
5; PURPOSE:
6;
7;       This function computes the autocorrelation Px(K,L) or
8;       autocovariance Rx(K,L) of a sample population X[nx,ny] as a
9;       function of the lag (K,L).
10;
11; CATEGORY:
12;       Statistics.
13;
14; CALLING SEQUENCE:
15;       Result = a_correlate2d(X, Lag)
16;
17; INPUTS:
18;       X:    an 2 dimension Array [nx,ny]
19;
20;     LAG:    2-element vector, in the intervals [-(nx-2), (nx-2)],[-(ny-2), (ny-2)],
21;             of type integer that specifies the absolute distance(s) between
22;             indexed elements of X.
23;
24; KEYWORD PARAMETERS:
25;       COVARIANCE:    If set to a non-zero value, the sample autocovariance
26;                      is computed.
27;
28;       DOUBLE:        If set to a non-zero value, computations are done in
29;                      double precision arithmetic.
30;
31; EXAMPLE:
32;
33; PROCEDURE:
34;
35;
36;                          nx-k-1  ny-l-1
37;                          sigma   sigma  (X[i,j]-Xmean)(X[i+k,j+l]-Ymean)
38;                           i=0     j=0
39;   correlation(X,[k,l])=------------------------------------------------------
40;                                  nx-1   ny-1                 
41;                                 sigma  sigma  (X[i,j]-Xmean)^2)
42;                                  i=0    j=0                     
43;
44;
45;                          nx-k-1  ny-l-1
46;                          sigma   sigma  (X[i,j]-Xmean)(Y[i+k,j+l]-Ymean)
47;                           i=0     j=0
48;   covariance(X,[k,l])=------------------------------------------------------
49;                                            nx*ny
50;
51;   Where Xmean is the mens of the sample population
52;   x=(x[0,0],x[1,0],...,x[nx-1,ny-1]).
53;
54;
55; REFERENCE:
56;
57; MODIFICATION HISTORY:
58;       28/2/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
59;       Based on the A_CORRELATE procedure of IDL
60;-
61
62FUNCTION Auto_Cov2d, X, Lag, Double = Double, zero2nan = zero2nan
63   XDim = SIZE(X, /dimensions)
64   nx = XDim[0]
65   ny = XDim[1]
66;Sample autocovariance function
67   Xmean = TOTAL(X, Double = Double) / (1.*nx*ny)
68;
69   res = TOTAL( (X[0:nx-1-lag[0], 0:ny-1-lag[1]] - Xmean) * $
70                (X[lag[0]:nx-1, lag[1]:ny-1] - Xmean) $
71                , Double = Double )
72   if keyword_set(zero2nan) AND res EQ 0 then res = !values.f_nan
73   RETURN, res
74
75END
76
77FUNCTION A_Correlate2d, X, Lag, Covariance = Covariance, Double = Double
78
79;Compute the sample-autocorrelation or autocovariance of (Xt, Xt+l)
80;as a function of the lag (l).
81
82   ON_ERROR, 2
83
84   XDim = SIZE(X, /dimensions)
85   XNDim = SIZE(X, /n_dimensions)
86   nx = XDim[0]
87   ny = XDim[1]
88   if XNDim NE 2 then $
89    MESSAGE, "X array must contain 2 dimensions."
90;Check length.
91   if nx lt 2 then $
92    MESSAGE, "first dimension of X array must contain 2 or more elements."
93   if ny lt 2 then $
94    MESSAGE, "second dimension of X array must contain 2 or more elements."
95   if n_elements(Lag) NE 2 THEN $
96    MESSAGE, "Lag array must contain 2 elements."
97   
98;If the DOUBLE keyword is not set then the internal precision and
99;result are identical to the type of input.
100   if N_ELEMENTS(Double) eq 0 then $
101    Double = (SIZE(X, /type) eq 5)
102
103   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Autocorrelation.
104      Auto = Auto_Cov2d(X, ABS(Lag), Double = Double) / $
105          Auto_Cov2d(X, [0L, 0L], Double = Double, /zero2nan)
106   endif else begin             ;Compute Autocovariance.
107      Auto = Auto_Cov2d(X, ABS(Lag), Double = Double) / n_elements(X)
108   endelse
109
110   if Double eq 0 then RETURN, FLOAT(Auto) else $
111    RETURN, Auto
112
113END
Note: See TracBrowser for help on using the repository browser.