3

What are some good ways of quantitatively characterizing the periodicity level of an approximately periodic $f(t), t \in \mathbb{R}$ signal?

I need this to tell if the output of some system is periodic or not, and to be able to decide the output of which system (out of a given set) is "more" periodic (approximates a perfectly periodic signal better).

I'm interested both in purely mathematical definitions for $-\infty < t < \infty$, and in practical ways of calculating this periodicity measure for a finite and time-discretized dataset. (Note that I expect that a useful measures might give higher periodicity for longer samples of the same signal for obvious reasons---but I'd like something that works for a series of any length.) The measure should (theoretically) be the same for infinitely long perfectly periodic signals, regardless of signal shape.

EDIT: Sample signals: http://ge.tt/91Zniy9

These signals are always between 0 and 1. I realized that the amplitude is also very relevant for my particular applications (very low amplitude signals should not be considered "periodic"), but any solution is easily modified to include this. I'm mentioning it to point out that solutions that are explicitly amplitude-dependent are also useful.

Szabolcs
  • 1,118
  • 1
  • 7
  • 27
  • Actually, _all_ examples are periodic ($f(t+T)=f(t)$). Periodicity has nothing to do with a function being jagged or not. Sines and cosines are not the only periodic functions! To distinguish between your functions you can try $\sum_t (f(t+1)-f(t))^2 / \langle (f-\langle f \rangle)^2 \rangle$ (it is normalized by the variance). For your particular examples it is 61.5, 62.0 and 1606.1. – Piotr Migdal Nov 15 '11 at 12:17

1 Answers1

2

Try using recurrence plot or perform a Fourier transform.

  • Recurrence plot

In the first case, the periodicity is translated to equally spaced lines with the slope 1 (and $R(t,t'+T) = R(t,t')$). This method is more robust, and easier to visualize, than the one below.

Then you can try to detect lines from the plot. E.g.:

  • Approach 1: Subtract the diagonal and calculate fractal dimension by box counting. The result, which is in the range of $[0,2]$, is 1 for a periodic signal, but not only for it.

  • Approach 2: Rotate by $45^\circ$ (or rather - take subimage rotated by $45^\circ$) and perform the Singluar Value Decomposition - then $\lambda_0/\sum_i \lambda_i$ (where $\lambda_0$ is the largest singular value) is an easy (but I guess not that robust) measure of periodicity (closer to one, the more periodic).

There are other methods for detection of lines (e.g. Radon transform) but here they seem to be using a canon to kill a fly.

  • Fourier transform

In the second case, for a periodic function its Fourier transform $\tilde{f}$ is a sum of equally spaced Dirac delta peaks: $$\tilde{f}(\omega) = \sum_{k=-\infty}^\infty a_k \delta(\omega - k \omega_0).$$ Of course for every practical data it will be smeared. If you want a quantification then $\int |\tilde{f}(\omega)|^4 d\omega$ may work (or use Shannon/Renyi entropy) - the larger it is, the more periodic is the function.

EDIT: Edited considerably to give numbers characterizing periodicity.

Piotr Migdal
  • 5,586
  • 2
  • 26
  • 70
  • Learning about recurrence plots was interesting, so I'm glad you mentioned them, but I'd like to point out that for this application what I need is a single real number to measure periodicity and make automatic decisions, not visualization. – Szabolcs Nov 11 '11 at 13:44
  • @Szabolcs: Then you can try to detect lines from the plot. E.g. 1) subtract the diagonal and calculate fractal dimension by box counting (the result should be 1 for periodic signal, but not only for it); 2) Rotate by 45deg (or rather - take subimage rotated by 45deg) and perform the [Singluar Value Decomposition](http://en.wikipedia.org/wiki/Singular_value_decomposition) - then $\lambda_0/\sum_i \lambda_i$ is an easy (but I guess not that robust) measure of periodicity (closer to one, the more periodic). – Piotr Migdal Nov 11 '11 at 13:56
  • Infinitely better answer than the standard "use Fourier transform" (-; –  Nov 11 '11 at 16:34
  • There are many possibilities starting from what you suggested, so this is a useful answer, thanks! However, the concrete suggestions are not the most practical for the problem at hand: I have a simulation written in C++ that after a while might start producing periodic output. I need to detect when the output is "periodic enough", then run the whole thing for $t$ more time to get enough data to analyse. While there are C++ libraries for what you suggested, the analysis phase that require libraries will be much more complex (and possibly slower) than the simulation itself, – Szabolcs Nov 12 '11 at 12:24
  • ... so I'd like to avoid this. Also, many external libraries will make it difficult for me to send the source code to others (read: my boss) and explain how to compile. So I'd prefer a simple and simple to implement solution (in C++, in Mathematica it's easy to implement). I'm looking at the recurrence plot now, perhaps I might be able to come up with something simple enough using that. – Szabolcs Nov 12 '11 at 12:27
  • @Szabolcs: Is it possible for you to post an example for a periodic signal and of an aperiodic signal? If you want to have a simple algorithm (i.e. without need to use any libraries) you need to specify very well how do you want to quantify aperiodicity. For example, for some applications $\sin(\omega(t) t)$ (with slowly changing period) is almost periodic (because locally it is), whereas for others - it is very aperiodic (as it has no global characteristic time). – Piotr Migdal Nov 12 '11 at 12:47
  • See my edit to the question, and sorry for the late reply! – Szabolcs Nov 15 '11 at 10:34
  • @PiotrMigdal Could you please elaborate (or provide some hints) on why you used a 4th power in your suggestion $\int |\tilde{f}(\omega)|^4 d\omega$? – Szabolcs Jan 16 '12 at 18:12
  • 1
    In short - $\int |f(t)|^2 dt = \int |\tilde{f}(\omega)|^2 d\omega$. While $f(t)$ and $\tilde{f}(\omega)$ are amplitudes, their abs val squares (intensities) can be treated as probabilities (given you normalize them, so they integrate to 1). Then one of way of saying if it is either 'peaked' or 'evenly distributed' is to calculate its entropy. One of variants is Renyi-Tsallis Entropy, which is in principle a function of $\int \int (|\tilde{f}(\omega)|^2)^\alpha$. I plugged $\alpha=2$ but you can play with other values. – Piotr Migdal Jan 16 '12 at 19:53