CODEWITHSUNDEEP
r
Emma Tebbs
Emma Tebbs

Reputation: 1467

estimate phase difference between two signals in r

If I have two time series, such as:

t <- seq(1,30)
y1 <- 4*sin(t*(2*pi/4) + 3)
y2 <- 4*cos(t*(2*pi/4) + 3)

plot(t,y1, type = 'l')
lines(t,y2, col = 'red')

I can calculate the phase difference between the signals in matlab using:

[Pxy,Freq] = cpsd(y1,y2);
coP = real(Pxy);
quadP = imag(Pxy);
phase = atan2(quadP,coP);

How would I achieve the same in R? What function can I use to estimate the phase difference between two time series in R

Upvotes: 2

Views: 3431

Answers (1)

RHertel
RHertel

Reputation: 23798

It is first necessary to determine the frequencies of the time series described by both datasets. If the frequencies turn out not to be equal, the notion of a phase difference doesn't make much sense. The frequency of the datasets can be extracted using the psd package.

#initial data
t <- seq(1,30)
y1 <- 4*sin(t*(2*pi/4) + 3)
y2 <- 4*cos(t*(2*pi/4) + 3)
# spectral analysis
library(psd)
out1 <- pspectrum(y1)
out2 <- pspectrum(y2)
f1 <- out1$freq[which.max(out1$spec)] # frequency with the highest peak in the spectrum
f2 <- out2$freq[which.max(out2$spec)]
# results:
#> f1
#[1] 0.25
#> f2
#[1] 0.25
f <- f1

This is a reassuring intermediate result. Firstly, the code has identified that for both time series the frequencies corresponding to the highest peak in the spectrum are equal. Secondly, the value of the frequency is now known, f=0.25. This is in agreement with the equations used to construct the dataset according to the OP, where a period of T=1/f=4 was chosen.

Both datasets, y1 and y2, can now be fitted to a function proportional to sin(2*pi*f*t)+cos(2*pi*f*t). The coefficients of these fits will provide information on the phase:

# fitting procedure:
fit1 <- lm(y1 ~ sin(2*pi*f*t)+cos(2*pi*f*t))
fit2 <- lm(y2 ~ sin(2*pi*f*t)+cos(2*pi*f*t))
#calculation of phase of y1:
a1 <- fit1$coefficients[2]
b1 <- fit1$coefficients[3]
ph1 <- atan(b1/a1)
#calculation of phase of y2:
fit2 <- lm(y2 ~ sin(2*pi*f*t)+cos(2*pi*f*t))
a2 <- fit2$coefficients[2]
b2 <- fit2$coefficients[3]
ph2 <- atan(b2/a2)
phase_difference <- as.numeric((ph2-ph1)/pi)
# result:
> phase_difference
#[1] 0.5

This means that the time series are phase-shifted by pi/2, as they should be according to the way the data has been generated.

For completeness, I include plots of the raw data and the fits:

> plot(y1~t)
> lines(fitted(fit1)~t,col=4,lty=2)

enter image description here

The blue and green dashed lines represent the function fitted to y1 and y2, respectively.

> plot(y2~t)
> lines(fitted(fit2)~t,col=3,lty=2)

enter image description here

Upvotes: 5

Related Questions