日长平滑变化

Question

我想对日长随时间平滑变化（但保持正弦曲线）的情况进行建模。用于更改瞬时频率的“啁啾”公式在https://en.wikipedia.org/wiki/Chirp 中给出，但在 5 天 24 小时内编码，然后再过 5 天过渡到 12 小时时，它看起来不正确：

period = list(  c(24,24,5), c(24,12,5) )
alpha = list(   c(0,5),     c(0,5)  )
s_samples = 100
A=50
O=50
simulatedData = data.frame(t=numeric(), v=numeric()) #initialise the output
daySteps = c(0, cumsum(unlist(period)[seq(3,length(unlist(period)), by=3)])) #set up the period starts and ends to set over, starting at 0
##Cycle over each of the items in the list
for(set in seq(period) ){
  t_points = s_samples*period[[set]][3]
  t = seq(daySteps[set], daySteps[set+1], length.out=t_points) #make the time
  slope = (24/period[[set]][2]-24/period[[set]][1])/(max(t)-min(t)) # get the slope
  f0 = 24/period[[set]][1] - slope*(min(t)) # find the freq when t0
  c = (24/period[[set]][2]-f0)/(max(t)) #calculate the chirp see https://en.wikipedia.org/wiki/Chirp and https://dsp.stackexchange.com/questions/57904/chirp-after-t-seconds
  wt = ((c*(t^2))/2) + f0*(t) # calc the freq 
  a = alpha[[set]][1]
  v = A * cos(2*pi*wt - a) + O
  simulatedData = rbind(simulatedData, data.frame(t, v) )
}
plot(simulatedData, type="l", lwd=2)
t = seq(0,sum(unlist(period)[seq(3,length(unlist(period)), by=3)]), by=1/24)
points(t, A*cos(2*pi*t)+O, col=3, type="l", lty=2)
points(t, A*cos(2*(24/12)*pi*t)+O, col=4, type="l", lty=2)

如预期的那样，前 24 天是完美的，第二个 5 天的最后一部分与 12 小时的循环相匹配，但该时期的第一部分看起来相差 180 度。怎么了？

Answer 1

我认为你让这变得比它需要的复杂得多。请记住，许多 R 函数已经向量化。以下函数将在t0 和t1 之间的频率f0 和f1 之间产生线性啁啾，并带有一个可选的phi 参数来指定您希望序列在周期的哪个点开始：

chirp <- function(f0, f1, t0, t1, phi = 0, n_steps = 1000)
{
  C <- (f1 - f0)/(t1 - t0)
  x <- seq(t0, t1, length.out = n_steps)
  y <- sin(2 * pi * (C / 2 * (x - t0)^2 + f0 * (x - t0)) + phi) # Ref Wikipedia
  data.frame(x, y)
}

当然，它也可以通过在两个相同频率之间“啁啾”来产生静态图的前半部分，因此我们可以通过做得到图上x，y点的数据框

df <- rbind(chirp(1, 1, 0, 5), chirp(1, 2, 5, 10))

结果：

plot(df$x, df$y, type = "l")

请注意，5 到 10 天之间有 7.5 个周期，因此如果您想顺利继续频率 2，则需要将 phi 参数设置为半周期（即 pi）：

df <- rbind(df, chirp(2, 2, 10, 15, phi = pi))

plot(df$x, df$y, type = "l")

请注意，如果啁啾信号发生在原始信号的偶数个周期内，啁啾信号和 2 Hz 信号的相位将仅在 n 秒后匹配。对于奇数，相位将偏离 180 度。这是线性啁啾的数学结果。为了看到这一点，让我们使用我们的函数来啁啾超过 6 秒，以便相位在 10 秒时匹配：

plot(df$x, df$y, type = "l")
lines(df2$x, df2$y, lty = 2, col = "green")
lines(df3$x, df3$y, lty = 2, col = "blue")
lines(df$x, df$y)