rtrigonometryfrequency-analysis

Smooth change of day length

I want to model what it might look like to have the day length change smoothly over time (but remain sinusoidal). The formula for a "chirp", to change the instantaneous frequency is given at https://en.wikipedia.org/wiki/Chirp but it doesn't look right when coded for a 24h period over 5 days and then a transition to 12h over another 5 days:

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)

black = simulated; blue = 12h period; green = 24h period

The first 24 are perfect, as expected, and the last part of the second 5 days matches a 12h cycled, but the first part of that period looks 180deg out of phase. What's wrong?

Solution

  • I think you're making this a lot more complex than it needs to be. Remember that many R functions are already vectorized. The following function will produce a linear chirp between frequencies f0 and f1 between t0 and t1, with an optional phi parameter to specify at what point on the cycle you want your sequence to begin:

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

    Of course, it can also produce the static first half of your plot by "chirping" between two identical frequencies, so we can get a data frame of x, y points on the plot by doing

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

    Which results in:

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

    enter image description here

    Note that between 5 and 10 days there are 7.5 cycles, so if you wanted to smoothly continue frequency 2, you would need to set the phi parameter to a half cycle (i.e. to pi):

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

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

    enter image description here

    Note that the phases of the chirped signal and a 2 Hz signal will only match after n seconds if the chirp occurs over an even number of periods of the original signal. For an odd number, the phase will be out by 180 degrees. This is a mathematical consequence of a linear chirp. To see this, let's use our function to chirp over 6 seconds so the phases match at 10 seconds:

    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)

    enter image description here