R help - Sep 2011 - Cross Spectrum : Conversion of 2-D spectrum into a single complex array

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')

Perhaps this question is inappropriate or posted to the wrong list?

Any guidance would be appreciated. Thanks.


--- On Fri, 9/23/11, Joseph Park <jpark.us at att.net> wrote:
> From: Joseph Park <jpark.us at att.net>
> Subject: Cross Spectrum : Conversion of 2-D spectrum into a single complex
array
> To: r-help at r-project.org
> Date: Friday, September 23, 2011, 5:02 PM
> 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')
> 
>