[R] estimate phase shift between two signals

[Dylan Beaudette]

On Wednesday 04 June 2008, Dylan Beaudette wrote:
> Hi,
>
> Are there any functions in R that could be used to estimate the phase-shift
> between two semi-sinusoidal vectors? Here is what I have tried so far,
> using the spectrum() function -- possibly incorrectly:
>
>
> # generate some fake data, normalized to unit circle
> x <- jitter(seq(-2*pi, 2*pi, by=0.1), amount=pi/8)
>
> # functions defining two out-of-phase phenomena
> f1 <- function(x) jitter(sin(x), amount=0.25)
> f2 <- function(x, a) jitter(sin(x + a), amount=0.25)
>
> # compute y-values
> # we are setting the phase shift arbitrarily
> s <- pi/1.5632198
> y1 <- f1(x)
> y2 <- f2(x, s)
>
>
> # plot:
> plot(x, y1, type='p', col='red', cex=0.5)
> lines(lowess(x, y1, f=0.25), col='red')
>
> points(x, y2, col='blue', cex=0.5)
> lines(lowess(x, y2, f=0.25), col='blue')
>
>
> # generate time series object
> comb.ts <- ts(matrix(c(y1, y2), ncol=2))
>
> # multivariate spectral decomposition
> spec <- spectrum(comb.ts, detrend=FALSE)
>
> # but how to interpret the phase estimate?
> mean(spec$phase)
>
> the mean 'phase' as returned from spectrum() does not seem to match the
> value used to generate the data... Am I mis-interpreting the use or output
> from spectrum() here? If so, is there a general procedure for estimating a
> phase-shift between two noisy signals? Would I first have to fit a smooth
> function in order to solve this analytically?
>
> Thanks in advance,

I should have thought about this a little bit more. Here is a brute-force 
method that may suffice for now, using nls() fit to sin(x + a).


# generate some fake data, normalized to unit circle
x <- jitter(seq(-2*pi, 2*pi, by=0.1), amount=pi/8)

# functions defining two out-of-phase phenomena
f2 <- function(x, a) jitter(sin(x + a), amount=0.25)

# compute y-values
# we are setting the phase shift arbitrarily
ps1 <- 1.5632198
ps2 <- 0.25

y1 <- f2(x, ps1)
y2 <- f2(x, ps2)


# plot:
plot(x, y1, type='p', col='red', cex=0.5)
lines(lowess(x, y1, f=0.25), col='red')

points(x, y2, col='blue', cex=0.5)
lines(lowess(x, y2, f=0.25), col='blue')


# generate time series object
comb.ts <- ts(matrix(c(y1, y2), ncol=2))

# multivariate spectral decomposition
spec <- spectrum(comb.ts, detrend=FALSE)

# but how to interpret the phase estimate?
mean(spec$phase)



# fit a simple sine function to each signal
fit1 <- nls(y1 ~ sin(x + a), start=list(a=0))
fit2 <- nls(y2 ~ sin(x + a), start=list(a=0))




# can we determine phase-shift by looking at zero-crossings?
# where function == 0 / changes sign

x.clean <- seq(-2*pi, 2*pi, by=0.01)


y1.clean <- predict(fit1, data.frame(x=x.clean))
y2.clean <- predict(fit2, data.frame(x=x.clean))


plot(x.clean, y1.clean, type='l', col='red')
lines(x.clean, y2.clean, type='l', col='blue')
points(x, y1, col='red', cex=0.5)
points(x, y2, col='blue', cex=0.5)
abline(h=0, lty=2)

#compute zero-crossings
y1.zero.idx <- which(abs(diff(sign(y1.clean))) == 2)
y2.zero.idx <- which(abs(diff(sign(y2.clean))) == 2)


points(x.clean[y1.zero.idx], y1.clean[y1.zero.idx], pch=16, col='red')
points(x.clean[y2.zero.idx], y2.clean[y2.zero.idx], pch=16, col='blue')


# how close are they?

# estimated
mean(x.clean[y2.zero.idx] - x.clean[y1.zero.idx])
[1] 1.3625


# original phase shift
(ps1 - ps2)
[1] 1.313220


the results appear to be good enough. Any thoughts on this approach, or ideas 
on a more elegant / proper implementation?

Cheers,



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341