Re: [R] estimate phase shift between two signals

From: ONKELINX, Thierry <Thierry.ONKELINX_at_inbo.be>
Date: Thu, 05 Jun 2008 11:56:21 +0200

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_at_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_at_r-project.org [mailto:r-help-bounces_at_r-project.org] Namens Dylan Beaudette
Verzonden: woensdag 4 juni 2008 20:23
Aan: r-help_at_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_at_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.

______________________________________________
R-help_at_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.
Received on Thu 05 Jun 2008 - 11:46:27 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 05 Jun 2008 - 12:30:36 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive