From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>

Date: Wed Aug 31 2011 - 13:07:35 MDT

Date: Wed Aug 31 2011 - 13:07:35 MDT

On 08/31/2011 12:48 PM, Anne Seidenglanz wrote:

*> Hi,
*

*>
*

*> I am new to NCL and used the 'specxy_anal' script to read in two time
*

*> series to calculate their cospectra as well as phase and coherence
*

*> spectra. In terms of analysis of the results, I was wondering if
*

*> anyone knew how exactly these spectra are calculated, especially the
*

*> phase and coherence, or where I can find information on this?
*

===================================================================

*>
*

*> In the script this is only stated as :
*

*>
*

*> 'spec = specxy_anal(time series 1,time series 2, d, sm, pct)'
*

*>
*

*> whereby 'spec_anal' calculates everything; phase, coherence etc. are
*

*> only attributes of 'spec' for which no further details are shown.
*

*>
*

*> d, sm, pct only refer to detrending, smoothing and tapering of the series
*

*>
*

=================================================================

d,sm,pct only refer to detrending, smoothing and tapering of the series

The underlying code is fortran 77.

The actual code is, well, not very clear.

(1) the cospectrum and quadrature spectrum are computed

(2) the output from (1) is smoothed

(3) the *smoothed* cospectrum and quadrature are used to compute the

coherence and phase

========================================================

c compute the cospectrum and quadrature spectrum

sclcs = sqrt( xinfo(6)*yinfo(6) )

cospc(1) = spcx(1)*spcy(1)

quspc(1) = spcx(1)*spcy(1)

do 60 n=0,nx2-1

cospc(n+2) = sclcs*

1 (work(kptrcx+n)*work(kptrcy+n)+work(kpticx+n)*work(kpticy+n))

60 quspc(n+2) = -sclcs* ! minus sign is to match nag output

1 (work(kpticx+n)*work(kptrcy+n)-work(kptrcx+n)*work(kpticy+n))

kopt = 1

javodd = 2*iabs(jave/2)+1

call swrnav (cospc(2),nspc-1,work(kptwgt)

1 ,javodd,kopt,spval,work,ier) ! can use work again

call swrnav (quspc(2),nspc-1,work(kptwgt)

1 ,javodd,kopt,spval,work,ier)

c compute the coherence and phase now that the co/quspc have been calculated

sclcoh = 1.0

do 70 n=1,nspc

phase(n) = atan2( quspc(n),cospc(n))*57.29578

if (spcx(n).ne.0. .and. spcy(n).ne.0.) then

coher(n) = sclcoh*(cospc(n)**2+quspc(n)**2)/(spcx(n)*spcy(n))

else

coher(n) = 0.

endif

70 continue

=========================================================

Some comments from the f77 code

c . output

c .

c . spcx,spcy - spectrum vector of length nspc for x and y series

c . spcx(1) - spectral estimate at freq = 0 (ignore)

c . spcx(nspc) - spectral estimate at freq = 0.5

c . cospc - cospectrum : vector of length nspc

c . this is the real part of the cross spectrum. it

c . measures the extent to which there are oscillations

c . with the same phase in the two series (or with opposite

c . sign ,i.e. , with a phase shift of half a cycle).

c . in other words : it measures the contribution of

c . different frequencies to the total cross-covariance

c . at zero lag.

c . quspc - quadrature spectrum : vector of length nspc

c . this is the imaginary part of the cross spectrum. it

c . measures the extent to which there are oscillations

c . with a phase difference of a quarter cycle in either

c . direction. i.e., it measures the contribution of

c . different frequencies to the total cross-covariance

c . of the series when all harmonics of one series

c . are delayed a quarter cycle relative to the other

c . relative to the other series.

c . coher - coherence squared : vector of length nspc

c . this is analogous to the square of the correlation

c . coef except that the coherence squared is a

c . function of frequency.

c . phase - phase (in degrees) : vector of length nspc.

c . a positive phase indicates that y leads x (i think!)

==============================================================

HTH ... Good Luck

_______________________________________________

ncl-talk mailing list

List instructions, subscriber options, unsubscribe:

http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Wed Aug 31 13:07:41 2011

*
This archive was generated by hypermail 2.1.8
: Wed Sep 07 2011 - 10:58:58 MDT
*