[R] Cross Spectrum : Conversion of 2-D spectrum into a single complex array
Joseph Park
jpark.us at att.net
Fri Sep 23 19:02:16 CEST 2011
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')
More information about the R-help
mailing list