Phase reference of a periodic signal Assume an arbitrary (discrete) signal that is periodic and known over a whole period. I need a way to select a characteristic point along the signal such that I can always retrieve it even when the signal is phase shifted.

Candidates could be zero crossings or extrema or inflections... provided they exist and are unique or can be discriminated from others.

For example, the position of the maximum would work well for a triangle wave, a little less accurately for sinusoid and would raise an ambiguity for $\sin x+\frac13\sin 3x$ as this one has two maxima of equal height.

The method should be robust to noise and work for relatively smooth signals.

Do you know any general feature/formula having the desired properties ?
Compute a DFT of your signal over a period of $T$ samples, extract the phase $\phi$ (in the $[0, 2\pi[$ range) of the first harmonic, and use the point $\frac{2\pi-\phi}{2 \pi}T$ as your landmark.

This will not correspond to a particular characteristic landmark in the time-domain waveform; but the procedure to retrieve it is by design invariant to phase shifts, distortion adding higher harmonics, and to a large extent noise.



If the waveform of your signal is known in advance you can also use it as a template, and select as a landmark the lag, in samples, at which the cross-correlation between your observed period and your template is maximal.

Python code here:

import numpy
import pylab

def shift(s, k):
  n = s.shape[0]
  return s[(numpy.arange(n) + k) % n]

T = 256
t = numpy.arange(0.0, T) / T

original = numpy.sin(t * 2 * numpy.pi) + numpy.sin(3 * t * 2 * numpy.pi) / 3
time_shifted = shift(original, 40)
distorted = shift(original + original ** 3, 200)
noisy = shift(original + numpy.random.randn(1, T).ravel() * 0.2, 180)

signals = [original, time_shifted, distorted, noisy]

for i, s in enumerate(signals):
  pylab.subplot(411 + i)
  pylab.plot(s)

  # Phase of first harmonic method.
  first_harmonic_coefficient = numpy.fft.rfft(s)[1]
  angle_2pi = -numpy.angle(first_harmonic_coefficient) + 2 * numpy.pi
  angle_2pi = angle_2pi % (2 * numpy.pi)
  landmark = numpy.round(angle_2pi / (2 * numpy.pi) * T)
  pylab.plot(landmark, 0, 'r+')

  # Cross-correlation method.
  xcorr = numpy.fft.irfft(numpy.fft.rfft(s) * numpy.fft.rfft(original))
  landmark = numpy.argmax(xcorr)
  pylab.plot(landmark, 0, 'g*')

pylab.savefig('plot.png')