Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

c_timecorrelate.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: 8.2 KB

Line
1	;+
2	; NAME:
3	; C_TIMECORRELATE
4	;
5	; PURPOSE:
6	; This function computes the "time cross correlation" Pxy(L) or
7	; the "time cross covariance" between 2 arrays (this is some
8	; kind of c_correlate but for multidimenstionals arrays) as a
9	; function of the lag (L).
10	;
11	; CATEGORY:
12	; Statistics.
13	;
14	; CALLING SEQUENCE:
15	; Result = c_timecorrelate(X, Y, Lag)
16	;
17	; INPUTS:
18	; X: an Array which last dimension is the time dimension of
19	; size n, float or double.
20	;
21	; Y: an Array which last dimension is the time dimension of
22	; size n, float or double.
23	;
24	; LAG: A scalar or n-element vector, in the interval [-(n-2), (n-2)],
25	; of type integer that specifies the absolute distance(s) between
26	; indexed elements of X.
27	;
28	; KEYWORD PARAMETERS:
29	; COVARIANCE: If set to a non-zero value, the sample cross
30	; covariance is computed.
31	;
32	; DOUBLE: If set to a non-zero value, computations are done in
33	; double precision arithmetic.
34	;
35	; EXAMPLE
36	; Define two n-element sample populations.
37	; x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
38	; y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95]
39	;
40	; Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7
41	; lag = [-5, 0, 1, 5, 6, 7]
42	; result = c_timecorrelate(x, y, lag)
43	;
44	; The result should be:
45	; [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
46	;
47	; PROCEDURE:
48	;
49	;
50	; FOR L>=0
51	;
52	; n-L-1
53	; sigma (X[k]-Xmean)(Y[k+L]-Ymean)
54	; k=0
55	; correlation(X,Y,L)=------------------------------------------------------
56	; n-1 n-1
57	; sqrt( (sigma (X[k]-Xmean)^2)*(sigma (Y[k]-Ymean)^2))
58	; k=0 k=0
59	;
60	;
61	;
62	; n-L-1
63	; sigma (X[k]-Xmean)(Y[k+L]-Ymean)
64	; k=0
65	; covariance(X,Y,L)=------------------------------------------------------
66	; n
67	;
68	; FOR L<0
69	;
70	;
71	; n-L-1
72	; sigma (X[k+L]-Xmean)(Y[k]-Ymean)
73	; k=0
74	; correlation(X,Y,L)=------------------------------------------------------
75	; n-1 n-1
76	; sqrt( (sigma (X[k]-Xmean)^2)*(sigma (Y[k]-Ymean)^2))
77	; k=0 k=0
78	;
79	;
80	;
81	; n-L-1
82	; sigma (X[k+L]-Xmean)(Y[k]-Ymean)
83	; k=0
84	; covariance(X,Y,L)=------------------------------------------------------
85	; n
86	;
87	; Where Xmean and Ymean are the time means of the sample populations
88	; x=(x[t=0],x[t=1],...,x[t=n-1]) and y=(y[t=0],y[t=1],...,y[t=n-1]),
89	; respectively.
90	;
91	;
92	;
93	; REFERENCE:
94	; INTRODUCTION TO STATISTICAL TIME SERIES
95	; Wayne A. Fuller
96	; ISBN 0-471-28715-6
97	;
98	; MODIFICATION HISTORY:
99	; - 01/03/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
100	; Based on the C_CORRELATE procedure of IDL
101	; - August 2003 Sebastien Masson
102	; update according to the update made in C_CORRELATE by
103	; W. Biagiotti and available in IDL 5.5
104	;-
105
106	FUNCTION TimeCross_Cov, Xd, Yd, M, nT, Ndim, Double = Double, ZERO2NAN = zero2nan
107	;Sample cross covariance function.
108
109	compile_opt hidden
110	;
111	case Ndim OF
112	1:res = TOTAL(Xd[0:nT - M - 1L] * Yd[M:nT - 1L] $
113	, Double = Double)
114	2:res = TOTAL(Xd[, 0:nT - M - 1L] Yd[*, M:nT - 1L] $
115	, Ndim, Double = Double)
116	3:res = TOTAL(Xd[, , 0:nT - M - 1L] * Yd[, , M:nT - 1L] $
117	, Ndim, Double = Double)
118	4:res = TOTAL(Xd[, , , 0:nT - M - 1L] Yd[, , *, M:nT - 1L] $
119	, Ndim, Double = Double)
120	ENDCASE
121	if keyword_set(zero2nan) then begin
122	zero = where(res EQ 0)
123	if zero[0] NE -1 then res[zero] = !values.f_nan
124	ENDIF
125	;
126	RETURN, res
127
128	END
129
130	FUNCTION C_Timecorrelate, X, Y, Lag, Covariance = Covariance, Double = Double
131
132	;Compute the sample cross correlation or cross covariance of
133	;(Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l).
134
135	ON_ERROR, 2
136
137	xsize = SIZE(X)
138	ysize = SIZE(Y)
139	nt = float(xsize[xsize[0]])
140	NDim = xsize[0]
141
142	if total(xsize[0:xsize[0]] NE ysize[0:ysize[0]]) NE 0 then $
143	MESSAGE, "X and Y arrays must have the same size and the same dimensions"
144
145	;Check length.
146	if nt lt 2 then $
147	MESSAGE, "Time dimension of X and Y arrays must contain 2 or more elements."
148
149	;If the DOUBLE keyword is not set then the internal precision and
150	;result are identical to the type of input.
151	if N_ELEMENTS(Double) eq 0 then $
152	Double = (Xsize[Xsize[0]+1] eq 5 or ysize[ysize[0]+1] eq 5)
153
154	if n_elements(lag) EQ 0 then lag = 0
155	nLag = N_ELEMENTS(Lag)
156
157	;Deviations
158	if double then one = 1.0d ELSE one = 1.0
159	Ndim = size(X, /n_dimensions)
160	Xd = TOTAL(X, Ndim, Double = Double) / nT
161	Xd = X - Xd[*]#replicate(one, nT)
162	Yd = TOTAL(Y, Ndim, Double = Double) / nT
163	Yd = Y - Yd[*]#replicate(one, nT)
164
165	if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
166
167	case NDim of
168	1:if Double eq 0 then Cross = FLTARR(nLag) else Cross = DBLARR(nLag)
169	2:if Double eq 0 then Cross = FLTARR(Xsize[1], nLag) else Cross = DBLARR(Xsize[1], nLag)
170	3:if Double eq 0 then Cross = FLTARR(Xsize[1], Xsize[2], nLag) $
171	else Cross = DBLARR(Xsize[1], Xsize[2], nLag)
172	4:if Double eq 0 then Cross = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $
173	else Cross = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag)
174	endcase
175
176	if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross Crossation.
177	for k = 0, nLag-1 do begin
178	if Lag[k] ge 0 then BEGIN
179	case NDim of
180	1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
181	2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
182	3: Cross[, , k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
183	4: Cross[, , *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
184	endcase
185	ENDIF else BEGIN
186	case NDim of
187	1: Cross[k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
188	2: Cross[*, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
189	3: Cross[, , k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
190	4: Cross[, , *, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
191	endcase
192	ENDELSE
193	ENDFOR
194	div = sqrt(TimeCross_Cov(Xd, Xd, 0L, nT, Ndim, Double = Double, /zero2nan) * $
195	TimeCross_Cov(Yd, Yd, 0L, nT, Ndim, Double = Double, /zero2nan))
196	Cross = temporary(Cross)/((temporary(div))[*]#replicate(one, nLag))
197	endif else begin ;Compute Cross Covariance.
198	for k = 0, nLag-1 do begin
199	if Lag[k] ge 0 then BEGIN
200	case NDim of
201	1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
202	2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
203	3: Cross[, , k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
204	4: Cross[, , *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
205	ENDCASE
206	ENDIF else BEGIN
207	case NDim of
208	1: Cross[k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
209	2: Cross[*, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
210	3: Cross[, , k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
211	4: Cross[, , *, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
212	ENDCASE
213	ENDELSE
214	endfor
215	endelse
216
217	if Double eq 0 then RETURN, FLOAT(Cross) else RETURN, Cross
218
219	END
220

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format