[R] Cross Spectrum : Conversion of 2-D spectrum into a single complex array

[Joseph Park]

Hi, I'm wondering why the spectrum() phase of quadrature 
couple isn't purely +/-pi. 

But mostly, I'm looking for a recommended way to take a 2-D
spectrum and convert it into a single complex array. 

Kindly consider:

# 10 Hz sine wave 10 seconds long sampled at 50 Hz
deltaT = 1/50
t      = seq(0, 10, deltaT)
w      = 2 * pi * 10

x = ts( sin( w * t ), deltat = deltaT )
y = ts( sin( w * t ), deltat = deltaT )

# Want the cross spectrum between x and y
Sxy = spectrum( cbind( x, y ), plot = F )

# The phase is perfectly zero
plot( Sxy $ freq, Sxy $ phase[ ,1], type = 'l' )

# Make y a quadrature component of x
y = ts( cos( w * t ), deltat = deltaT )
Sxy = spectrum( cbind( x, y ), plot = F )

# The phase should be either +pi or -pi
# since exp(i pi) = exp(-i pi) = -1 + 0i
# Why isn't it?  Sampling issues?  Nyquist has been satisfied.
plot( Sxy $ freq, Sxy $ phase[ ,1], type = 'l' )

# The real question (limited to a 2-D input):
# How to combine the amplitude/phase into a single
# complex valued cross spectrum?
mySxy = complex( real     = Sxy $ spec[,1] * Sxy$ spec[,2],
                imaginary = Sxy $ phase[ ,1] )

# This does give something affine to a plot of Sxy
# Compare
plot( Sxy, log="dB" )
# to:
plot( Sxy$freq, 10*log10( Re(mySxy) ), type='l')