r - 估计 r 中两个信号之间的相位差

Question

如果我有两个时间序列，例如:

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')

我可以使用 matlab 计算信号之间的相位差:

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

我如何在 R 中实现同样的目标？我可以使用什么函数来估计R中两个时间序列之间的相位差

Answer 1

首先需要确定两个数据集描述的时间序列的频率。如果频率不相等，则相位差的概念就没有多大意义。可以使用 psd 包提取数据集的频率。

#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

这是一个令人放心的中间结果。首先，代码确定对于两个时间序列，频谱中最高峰值对应的频率是相等的。其次，现在已知频率值，f=0.25。这与根据 OP 构建数据集所使用的方程一致，其中选择的周期为 T=1/f=4。

两个数据集 y1 和 y2 现在都可以拟合到与 sin(2*pi*f*t)+cos(2 *pi*f*t)。这些拟合的系数将提供有关相位的信息:

# 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

这意味着时间序列相移 pi/2，因为它们应该根据数据的生成方式进行相移。

为了完整起见，我添加了原始数据和拟合的图表:

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

蓝色和绿色虚线分别表示拟合到 y1 和 y2 的函数。

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