估计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)