[R] estimate phase shift between two signals

[ONKELINX, Thierry]

Dylan,

You can write r sin(x + alfa) as a sin(x) + b cos(x). Then you just have
to fit a linear model. r =  sqrt(a^2 + b^2) and tan(alfa) = b / a

HTH,

Thierry


------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
Thierry.Onkelinx op inbo.be 
www.inbo.be 

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-----Oorspronkelijk bericht-----
Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org]
Namens Dylan Beaudette
Verzonden: woensdag 4 juni 2008 20:23
Aan: r-help op r-project.org
Onderwerp: Re: [R] estimate phase shift between two signals

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

______________________________________________
R-help op r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.