Re: specxy_anal

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
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