Finding the local characteristic period of a signal

Question

Sorry if this is a basic question, but suppose I have a signal that is locally periodic over short time intervals but is not globally of that period. As an example:

import numpy as np
import matplotlib.pyplot as plt

x_limit = 100
res = 50000

random_factor = 0.345
periods_per_slip = 4

x = np.linspace(-x_limit, x_limit, res)
y = np.sin(2*np.pi*x)
yp = 1.1*np.sin(2*np.pi*x - random_factor * (np.pi * np.floor(x/periods_per_slip)))

plt.plot(x,y)
plt.plot(x, yp)
plt.xlim(-10, 10)

The blue curve is a pure sinusoid of period 1, while the orange curve is locally a sinusoid of the same period but with occasional phase slips.

The Fourier transform is below

plt.figure()
plt.plot(np.linspace(-np.pi/(2*x_limit/res), np.pi/(2*x_limit/res), res), np.abs(np.fft.fftshift(np.fft.fft(y))))
plt.plot(np.linspace(-np.pi/(2*x_limit/res), np.pi/(2*x_limit/res), res), np.abs(np.fft.fftshift(np.fft.fft(yp))))
plt.xlim(-10, 10)
plt.ylim(-1000, 30000)

Here, the dominant frequency is not $\omega = 2\pi$ (as it is in the pure sinusoid) but is shifted slightly.

So my question is: what is a proper analysis method for extracting the local characteristic period ($2\pi/\omega = 1$) of this signal? Short-time Fourier transforms? Wavelets? Clearly, just taking the frequency with the highest FFT amplitude gives the "wrong" answer.

The example above is just a toy model. My actual application is image processing for crystallography, so I'm looking for a technique that generalizes to two dimensions.

Answer 1

Clearly, just taking the frequency with the highest FFT amplitude gives the "wrong" answer.

That's a matter of interpretation. The main reason why the peak has shifted is that all your phase slips are a going in the same direction. Hence your average rate of phase accumulation is lower and hence the dominant frequency is lower. If you had chosen mean free phase slips, the main frequency would have been preserved.

what is a proper analysis method for extracting the local characteristic period

That depends A LOT on how exactly you define the characteristic period and what is the mechanism of the period change is.

You can always try chopping up the signals into frames that have a length that's similar to the time scale or your expected changes. Than you can try short term Fourier Transform or running cross correlation. Another good option would be a PLL. You still would have to define how you handle a frame with a single phase slip in it. What exactly IS the correct answer in this case?

If you know, that you are getting phase slips, you can write and algorithm that detects and discards them. Again, a PLL would be useful for this as the phase slip would generate a big jump in the error signal.

Answer 2

A real alternative (in some sense) to the Fourier Transform is the Hilbert Transform (HT), which enables amplitude-phase demodulation of your signal and computation of the instantaneous frequency by differenciation of the phase. In Matlab, your even don't need to perform the HT: just call instfreq on your signal.

Caution about the instantaneous frequency: its definition and existence has been a subject of much debate. I can give you references if needed.