The Nyquist-Shannon sampling theorem tells us the frequency above which periodic signals are not detectable in regularly-sampled timeseries measurements. For example, if we observe a star at regular intervals of $\Delta t = 1 \rm{hr}$, then the sampling frequency is $f_s = 1 / \Delta t$ and the Nyquist-Shannon sampling theorem tells us that the data stream is sensitive only to frequencies below the band-limit $B = f_s / 2$. Stated in another way: if our observing cadence is $\Delta t$, the shortest period we can hope to detect within the data is $2\Delta t$.
The above holds for regularly-sampled data. What about irregularly-sampled data? A common assumption seems to be that the operative criterion in the above proof is the minimum time $\Delta t$ between samples, so that for irregularly-sampled data at times $t_i$ the band-limit is at $B = 2\min_i(t_{i + 1} - t_{i})$. But this is not the case! What is more important is the window function of observations; as we'll see, when the sampling is irregular, it is possible to detect periods far smaller than the minimum $\Delta t$.
To understand why Nyquist doesn't apply in irregular data, we must first understand Fourier transforms and how they relate to this problem. The Fourier transform can be thought of, roughly, as a transormation expressing any function $h(t)$ into a frequency-space representation $H(f)$. This mapping is one-to-one and reversible, so that $H(f)$ and $h(t)$ in some senses contain the exact same information. You can read about the math elsewhere, but there are a few important features of the transform that it helps to understand conceptually.
This is the first important point to understand. Fourier analysis is constructed around sinusoids as the fundamental building block or basis function of periodic signals. Because of this, a simple sinusoidal signal maps to a delta-function in frequency space:
In [311]:
plot_sinusoid_FT() # function defined below
In [313]:
plot_sinusoid_sum_FT()
In [312]:
plot_nonperiodic()
In [342]:
plot_nonperiodic_2()
In [314]:
plot_windowed_sine()
In [327]:
plot_windowed_sine_sum()
In [252]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# Seaborn sets some nice plotting defaults.
# This can be safely commented-out
import seaborn; seaborn.set()
# AstroML contains some tools for approximating
# Continuous Fourier transforms. We'll use it below.
from astroML import fourier
In [253]:
def plot_sinusoid_FT():
"""Plot a single sinusoid and its Fourier power"""
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
func = lambda t: np.sin(2 * np.pi * t)
t = np.linspace(0, 5, 1000)
ax[0].plot(t, func(t))
ax[0].set(ylim=(-1.2, 1.2), xlabel='$t$', ylabel='$h(t)$',
title='Input Data')
t = np.linspace(0, 20, 10000)
f, H = fourier.FT_continuous(t, func(t))
ax[1].plot(f, abs(H))
ax[1].set(xlim=(-3, 3), xlabel='$f$', ylabel='$|H(f)|$',
title='Spectral Power')
plot_sinusoid_FT()
In [328]:
def plot_sinusoid_sum_FT():
"""Plot a sum of sinusoids and their Fourier power"""
def func(t, coeffs):
frequencies = 2 * np.pi * np.arange(1, len(coeffs) + 1)
return np.dot(coeffs, np.sin(np.outer(frequencies, t)))
coeffs = [1, 0.5, 0.25]
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
t = np.linspace(0, 3, 1000)
ax[0].plot(t, func(t, coeffs))
ax[0].set(ylim=(-1.8, 1.8), xlabel='$t$', ylabel='$h(t)$',
title='Input Data (sum of 3 sinusoids)')
for i, coeff in enumerate(coeffs):
ax[0].plot(t, coeffs[i] * np.sin(2 * np.pi * (i + 1) * t),
color='gray', alpha=0.3)
t = np.linspace(0, 20, 10000)
f, H = fourier.FT_continuous(t, func(t, coeffs))
ax[1].plot(f, abs(H))
ax[1].set(xlim=(-5, 5), xlabel='$f$', ylabel='$|H(f)|$',
title='Spectral Power')
plot_sinusoid_sum_FT()
In [306]:
def plot_nonperiodic():
window = lambda t: (abs(t - 2.5) < 1.5).astype(float)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
t = np.linspace(0, 5, 1000)
ax[0].plot(t, window(t))
ax[0].set(ylim=(-0.2, 1.5),
xlabel='$t$', ylabel='$h(t)$',
title='Input Data')
t = np.linspace(-50, 50, 10000)
f, H = fourier.FT_continuous(t, window(t))
ax[1].plot(f, abs(H))
ax[1].set(xlim=(-5, 5), xlabel='$f$', ylabel='$|H(f)|$',
title='Spectral Power')
plot_nonperiodic()
In [341]:
def plot_nonperiodic_2():
window = lambda t: (abs(t - 2.5) < 1.5).astype(float)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
t = np.linspace(-50, 50, 10000)
f, H = fourier.FT_continuous(t, t)
H = np.exp(-abs(abs(f) - 3) / 0.3)
t, h = fourier.IFT_continuous(f, H)
ax[0].plot(t, h)
ax[0].set(xlim=(-5, 5), #ylim=(-0.2, 1.5),
xlabel='$t$', ylabel='$h(t)$',
title='Input Data')
ax[1].plot(f, abs(H))
ax[1].set(xlim=(-10, 10), ylim=(0, 1.2), xlabel='$f$', ylabel='$|H(f)|$',
title='Spectral Power')
plot_nonperiodic_2()
In [310]:
def invisible_ticks(ax):
for label in ax.get_xticklabels():
label.set_visible(False)
def plot_windowed_sine():
window = lambda t: (abs(t - 2.5) < 1.5).astype(float)
func = lambda t: np.sin(2 * np.pi * 2 * t)
fig = plt.figure(figsize=(12, 8))
gs = plt.GridSpec(4, 2)
# plot inputs
t = np.linspace(0, 5, 1000)
ax = fig.add_subplot(gs[0, 0])
ax.plot(t, func(t))
ax.set(ylim=(-1.2, 1.2))
invisible_ticks(ax)
ax.set_title('Multiplication in time-space')
ax = fig.add_subplot(gs[1, 0], sharex=ax)
ax.plot(t, window(t))
ax.set(ylim=(-0.2, 1.2))
invisible_ticks(ax)
ax = fig.add_subplot(gs[2:, 0], sharex=ax)
ax.plot(t, func(t) * window(t))
ax.set(ylim=(-1.2, 1.2),
xlabel='$t$', ylabel='$h(t)$')
t = np.linspace(-50, 50, 10000)
ax = fig.add_subplot(gs[0, 1])
f, H = fourier.FT_continuous(t, func(t))
ax.plot(f, abs(H))
ax.set_ylim(0, 60)
invisible_ticks(ax)
ax.set_title('Convolution in frequency-space')
ax = fig.add_subplot(gs[1, 1], sharex=ax)
f, H = fourier.FT_continuous(t, window(t))
ax.plot(f, abs(H))
invisible_ticks(ax)
ax = fig.add_subplot(gs[2:, 1], sharex=ax)
f, H = fourier.FT_continuous(t, func(t) * window(t))
ax.plot(f, abs(H))
ax.set(xlim=(-6, 6), xlabel='$f$', ylabel='$|H(f)|$')
plot_windowed_sine()
In [326]:
def invisible_ticks(ax):
for label in ax.get_xticklabels():
label.set_visible(False)
def plot_windowed_sine_sum():
window = lambda t: (abs(t - 2.5) < 1.5).astype(float)
def func(t, coeffs=(1, 0.5, 0.25)):
frequencies = 2 * np.pi * np.arange(1, len(coeffs) + 1)
return np.dot(coeffs, np.sin(np.outer(frequencies, t)))
fig = plt.figure(figsize=(12, 8))
gs = plt.GridSpec(4, 2)
# plot inputs
t = np.linspace(0, 5, 1000)
ax = fig.add_subplot(gs[0, 0])
ax.plot(t, func(t))
ax.set(ylim=(-1.5, 1.5))
invisible_ticks(ax)
ax.set_title('Multiplication in time-space')
ax = fig.add_subplot(gs[1, 0], sharex=ax)
ax.plot(t, window(t))
ax.set(ylim=(-0.2, 1.2))
invisible_ticks(ax)
ax = fig.add_subplot(gs[2:, 0], sharex=ax)
ax.plot(t, func(t) * window(t))
ax.set(ylim=(-1.5, 1.5),
xlabel='$t$', ylabel='$h(t)$')
t = np.linspace(-50, 50, 10000)
ax = fig.add_subplot(gs[0, 1])
f, H = fourier.FT_continuous(t, func(t))
ax.plot(f, abs(H))
ax.set_ylim(0, 60)
invisible_ticks(ax)
ax.set_title('Convolution in frequency-space')
ax = fig.add_subplot(gs[1, 1], sharex=ax)
f, H = fourier.FT_continuous(t, window(t))
ax.plot(f, abs(H))
invisible_ticks(ax)
ax = fig.add_subplot(gs[2:, 1], sharex=ax)
f, H = fourier.FT_continuous(t, func(t) * window(t))
ax.plot(f, abs(H))
ax.set(xlim=(-6, 6), xlabel='$f$', ylabel='$|H(f)|$')
plot_windowed_sine_sum()
In [319]:
def plot_regular_window():
window = lambda t: (abs(t - np.round(t)) < 0.01).astype(float)
func = lambda t: np.sin(2 * np.pi * t)
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
t = np.linspace(0, 5, 1000)
ax[0].plot(t, func(t) * window(t))
ax[0].set(ylim=(-1.5, 1.5),
xlabel='$t$', ylabel='$h(t)$',
title='Input Data')
t = np.linspace(-50, 50, 10000)
f, H = fourier.FT_continuous(t, func(t) * window(t))
ax[1].plot(f, abs(H))
ax[1].set(xlim=(-6, 6), xlabel='$f$', ylabel='$|H(f)|$',
title='Spectral Power')
plot_regular_window()
In [ ]: