Estimating time periods of oscillators

[Sh_166]

As the title says, I want to estimate the time periods of a certain kind of oscillation that comes out of an ODE solution:

This behaviour has a large period among other smaller ones (in this image, the time period of the slow modulation is around 70 units) and I want to estimate this time period numerically.

I cannot use a visual estimate (like I did here) because I want to vary a certain parameter (F) in the problem and find out how the slowly varying period depends on the said parameter - so, I need to solve the DE in a loop. If I had a way to estimate the time period from the ODE solution set, I could’ve saved it into a vector or array for each run of the loop and plotted them afterwards to get the dependence.

I tried estimate_period() from ChaosTools.jl with lombscargle initially - but this gives me the smaller periods and not the one I want. I know Fourier transforms are usually used to decompose signals but I am not in the field of signal processing and am unfamiliar with this kind of usage.

I tried the FFTW.jl package but its going over my head (mostly because the documentation seems to lack examples) - even then I tried a simple code:

using Symbolics, Plots, DifferentialEquations, ModelingToolkit, LaTeXStrings, Measures, EasyModelAnalysis, CurveFit, ChaosTools, DSP, FFTW

@variables t x(t) y(t)
@parameters F, ω
d = Differential(t)

E₁ = d(d(x)) + ω^2*sin(x)
E₁₁ = Symbolics.solve_for(E₁ ~ ω^2*F*cos(ω*t), d(d(x)))
system = [d(d(x)) ~ E₁₁]

@named model = ODESystem(system, t, [x], [F, ω])

model = structural_simplify(model)

X = ODEProblem(model, u_1, (0, 100), [ω => 1, F => 0.1])
Y = solve(X, Rodas4P())
v1 = Y[2,:]
v2 = Y[1,:]

Fal = fft(v1) |> fftshift
freqs = fftfreq(length(Y.t)) |> fftshift
plot(freqs, @. abs(Fal))

This gives me:

I have no idea which of these peaks I want (I’d guess the ones closest to 0 since that should have the largest time period - but that is around 0.003 which gives me a time period of 1/0.003 \sim 333 units which is completely wrong). The farthest peaks are correspond to 100 time units as well - much larger than the visual estimate. So, my code must be wrong.

I would like to know if there’s a way to systematically do something like this - any help would be really appreciated!

[PetrKryslUCSD]

I think you have beats caused by the interaction of two signals of similar frequencies. The beat “frequency” is the difference.

[Benny]

The difficulty here is that the modulation frequency won’t directly show up as peak frequency but the difference between some peak frequencies, and the number of peaks is affected by whether a constant was added to the modulation sinusoid. Eyeballing the jagged graph again, I’m guessing 2 peaks at 0.0045 and 0.006, the difference being 0.0015 with a reciprocal of 667, again much larger than the visual estimate. Is fftfreq off? It defaults to a sampling rate of 1 sample per time unit, which seems lower than that from eyeballing the x(t) graph. A lower sampling rate than actual would stretch to longer periods, thus lower frequencies.

Incidentally, that link’s question also highlights that it’s important not to just eyeball the modulated waveform’s width as a simple sinusoid because the width involves an absolute value. Adding a constant 1 made the modulation sinusoid positive, so an absolute value does not affect it, and the width peaks once per cycle. However, an absolute value of a modulation sinusoid centered on 0 turns the valleys to peaks, so the width peaks twice per cycle.

[fatteneder]

Side note: fft assumes that your time series is sampled with a constant rate. You should check if that is true for the ODE solution returned from solve, but I suspect that it uses some adaptive time stepping.
To fix that you can interpolate your solution to a equidistant grid before applying the fft, e.g.

t = range(0, Y.t[end], length=1000)
YY = stack(Y.(t),dims=2)
v1 = YY[1,:]
v2 = YY[2,:]

Alternatively, you could instruct the ODE solver to not use adaptive stepping and give it a fixed step width, e.g. (untested)

Y = solve(X, Rodas4P(), adaptive=false, dt=1)

[Benny]

Oh since you’re trying out FFTs, there are a few general tips, especially since your FFT turned out a tad jagged. At this point, you could

1. zero-pad your input (just extending the vector with trailing 0s) to increase the frequency resolution of the FFT, which may let you estimate the peak frequencies better. This doesn’t actually add information, it’s rooted in how the discrete Fourier transform and its input is inherently periodic, whereas your x(t) is likely a nonperiodic window. Zero-padding will approach that.

2. If you zero-pad, you’ll notice more peaks trailing off the big ones in the magnitude spectrum (abs). This involves the sinc function, which is what happens when you take the continuous Fourier transform of a rectangular window, which is what zero-padding is approaching. For a continuous Fourier transform, a nice window is a smooth one-peak Gaussian function because its transform is another Gaussian, but the infinitely long domain can’t be represented digitally, and a well-approximated truncation is often undesirably narrow. You don’t want a narrow window because it erases information, widening the frequency peaks, but all windows will make a tradeoff between primary peak width and trailing peak height, and you’ll have to decide which you want. Note that applying the window before zero-padding is proper. DSP.jl provides some common windows, though not the Blackman-Harris window which I recall to be a common personal choice for short trailing peaks with a moderately widened primary peak.

If you want to dive further into why these tips exist, I tracked down this stackexchange post with a lot of math relating the continuous and discrete Fourier transforms and explaining the general tips. It doesn’t quite spell out how the continuous transform of the Dirac-sampled and truncated signal is related to the continuous transform of the original signal when you don’t sample full periods (like this case), and I’ll need to wrap my head around that a little.

So the first derivation block’s 1st line shows the continuous transform of the original signal being convolved with a sinc function of height 1 and repeated over its cyclic frequency domain; this is equivalent to the 3rd line expressing a 1/M-scaled continuous transform of the rectangular windowed and Dirac-sampled signal. The preceding one-sentence description does not involve the 1st line or the 1/M scaling. As a later bullet point mentions, the 1/M scaling normalizes the magnitude of the 3rd and last lines to that of the 1st line. What isn’t spelled out is that the 1st line’s magnitude is good to normalize to because it does not vary over M, and the sinc peak heights can then estimate the area under any Dirac delta peaks in the original signal’s magnitude spectrum (from sinusoids); you adjust the estimate for different windows. This isn’t a widely applicable strategy because it only works well for Dirac delta peaks, not other curves.

It may seem simpler to deconvolve the 1st line to retrieve the original signal’s transform, but that’s technically not possible because it would be equivalent to dividing by 0 outside the window and magically restoring those parts of the signal. For example, consider a gaussian-windowed sinusoid, a kind of wavelet, that has a transform with gaussian peaks at the sinusoidal frequency. If that wavelet is outside the window, that gaussian frequency peak is erased. So, however you try to estimate the original signal, including some deconvolution approaches, you need to be aware that the sampled window is limited information. If the window is a segment of a periodic or otherwise unchanging wave, it could be a good estimate. I’m not well-versed enough to offer more, but there are packages Deconvolution.jl and Peaks.jl pertaining to what was said here.

[mbaz]

It looks like the signal is the product of two sinusoids:

\cos(2\pi f_\text{small}t) \cos(2\pi f_\text{large} t).

(Omitting the angle for clarity). This is equal to

\cos(2\pi (f_\text{large} - f_\text{small}t) + \cos(2\pi (f_\text{large} + f_\text{small}t).

There are some tools that can be used to estimate those two frequencies:

• DSP ESPRIT implementation
• There are a few other frequency estimators in this PR.
• There are many related questions/answers in DSP.stackexchange; one example.

[Sh_166]

Thanks a lot to everyone in the thread for their inputs - I know nearly nothing about signal analysis, so it helps a lot! I haven’t yet gone through the links posted by @Benny - but they might help in understanding the issue.

I think my mistake was the sampling rate issue - the ODE solver uses an adaptive timestep. I used @fatteneder’s tip to sample it evenly and recreate the vectors v1, v2 and used the ESPRIT implementation mentioned by @mbaz. This gave me two frequencies - one just the negative of the other and with magnitude equal to 0.0149 which gives me roughly a period of 68 time units.

The modified MWE is:

using Symbolics, Plots, DifferentialEquations, ModelingToolkit, LaTeXStrings, Measures, EasyModelAnalysis, CurveFit, ChaosTools, DSP, FFTW

#Same as last MWE

X = ODEProblem(model, u_1, (0, 100), [ω => 1, F => 0.1])
Y = solve(X, Rodas4P())

t = range(0, Y.t[end], length=1000)
YY = stack(Y.(t),dims=2)
v1 = YY[1,:]
v2 = YY[2,:]

esprit(v1, 333, 2) #length, N = 1000, so M = (N+1)/3 ~ 333

I am not sure I understand the code very well - for a start, the ESPRIT documentation mentions that the default sampling rate for the function is 1.0 Hz while my t span is (0, 100) (so, a length of 1000 leads to a 0.1 interval - or, a 0.1 time units sampling rate I think which is different to the documentation). Changing the length of the t vector to 100 (and the M value to 33) gives me a period of 6.7 time units which is too small.

However, this is also more or less the period I expect for the faster oscillation - so, I am not sure if that is what I found here.

[rafael.guerra]

It may be easier to find the lower frequency if you run the simulation for at least one period (140 s).

See the proxy code below if it helps.

Proxy example

using Plots, FFTW
Δt = 0.1            # sampling period [s]
t = 0:Δt:140
ω₀ = 2π/140         # angular frequency [Hz]

A(t) = cos(20*ω₀*t - π/2)
m(t) = A(t)*cos(ω₀*t + π/2)

yt = m.(t)
p1 = plot(t, yt)

Yf = fft(yt)
f = fftfreq(length(yt), 1/Δt)   # frequencies of Yf

# Keep positive frequencies:
n = length(t)
ns = isodd(n) ? (n-1)÷2 : n÷2-1 
AYf = abs.(Yf)[1:ns]
f = f[1:ns]

p2 = plot(f, AYf, xlims=(0, 0.3))
plot(p1, p2)


ix = partialsortperm(AYf, 1:2, rev=true)
f1, f2 = f[ix]
Δf = abs(f2 - f1)       # 0.0143 Hz
T₀ = 2/Δf               # 140 s

[Benny]

@Sh_166

I am not sure I understand the code very well - for a start, the ESPRIT documentation

Another aspect is that ESPRIT requires you know your signal is a sum of sinusoids (in your code, 2), but as the 1st stackexchange post I link showed, a modulated sinusoid can be >2. Plotting the magnitude of a DFT after windowing and zero-padding can help confirm the appropriate number. For frequencies so close to each other, you might have to try out several windows because wider main peaks and taller trailing peaks can overlap and distort the peak maxima away from the frequencies, but you only need it to be clear enough to identify the number of sinusoids for ESPRIT.

Note that these difficulties with close peaks aren’t sampling limitations, but trigonometry. For an identity pertaining to modulated sinusoids, cos(m)cos(f) = [cos(f-m) + cos(f+m)]/2, and we accordingly perceive 2 pure tones with close pitches as 1 intermediate pitch that periodically varies in loudness.

[ChrisRackauckas]

@fatteneder

Alternatively, you could instruct the ODE solver to not use adaptive stepping and give it a fixed step width, e.g. (untested)

No, please don’t change the ODE solver stepping in order to change the output. This is turning off error estimation so there’s no guarantees on convergence. Use saveat or post-solution interpolation. Not only is it simpler, saveat=1 will give you evenly spaced output for the FFT, it’s also more numerically correct.

[StackEnjoyer]

Maybe this might help

"""
    find_periods(data, sw::Int=2000, step::Int=sw÷2)

Finds period indices (minima) within a time series `data` using a sliding window approach.
The method is most robust when the minima are distinguishable from noise and
when the window width `sw` is shorter than the expected period in the data.

# Arguments
- `data`: The input time series data.
- `sw`: The width of the sliding window (default: 2000).
- `step`: The step size for the sliding window (default: sw÷2).

# Returns
- `period_indices`: An array of period indices detected within the input data.
"""
function find_periods(data, sw::Int=2000, step::Int=sw ÷ 2)
    period_indices = Int64[]  # Initialize an array to store the detected period indices
    min_counter = 0           # Initialize a counter for consecutive minima
    min_index = 1             # Initialize the index of the first minimum
    trigger = sw ÷ step - 1     # Define the trigger value based on the window width and step size
    # Iterate through the data using the sliding window approach
    for i = 1:step:length(data)-sw
        buffer = @view data[i:i+sw]
        ~, current_index = findnanmin(buffer)
        # Check if the current minimum is the same as the previous minimum
        if min_index == current_index + i
            # Increment the counter for consecutive minima
            min_counter += 1
        else
            # Update the minimum index and reset the counter
            min_index = current_index + i
            min_counter = 0
        end
        # Is a period minima if the counter reaches the trigger value
        if min_counter == trigger
            push!(period_indices, min_index)
        end
    end
    return period_indices
end

[Sh_166]

Thanks to everyone for their help - much appreciated! I figured out a simpler but slightly more inaccurate way of doing this, but this is dependent completely on my problem.

I was already using EasyModelAnalysis.jl to find the maximum amplitude of the oscillation via get_max_t() and this function returns both the maximum amplitude and the time point it occurs at. I figured that simply multiplying this time by 2 would give me a very simple (albeit not very accurate all the time) estimate of the time period I wanted to extract. This is exactly what I did and it more or less worked.

Thanks again to everyone for their help!

[rafael.guerra]

@Sh_166

get_max_t() and this function returns both the maximum amplitude and the time point it occurs at. I figured that simply multiplying this time by 222 would give me a very simple (albeit not very accurate all the time) estimate of the time period

In the OP, this gives a time around 35 s. You need to multiply by 4 to get the (approximate) period of the low frequency carrier, not by 2.

[Sh_166]

I meant the time for one pulse - that extends from 0 to 70 s approximately.