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 [EMAIL PROTECTED] 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens Dylan Beaudette Verzonden: woensdag 4 juni 2008 20:23 Aan: r-help@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@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@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.