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

[Joseph Park]

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