Phase reference of a periodic signal
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')