日长的平滑变化

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)

日长的平滑变化

Smooth change of day length

trigonometry

r

frequency-analysis