Usage#
fourierror propagates measurement uncertainty through a discrete Fourier transform. It offers four entry points, and the point of this notebook is to show what each one is for and where they agree:
Function |
Method |
Cost |
|---|---|---|
|
exact propagation through the DFT matrix |
O(N³) matrix products |
|
the same exact result, computed with FFTs |
O(N² log N) |
|
Monte Carlo sampling of the input distribution |
|
|
Monte Carlo sampling, via the FFT |
|
Everything is built on scipp data objects, so inputs carry units, coordinates and variances, and the frequency axis is derived for you.
[1]:
import numpy as np
import scipp as sc
import matplotlib.pyplot as plt
from fourierror import analytical, numerical
from fourierror.variable import CovVariable
from fourierror.data_array import CovDataArray
Input data#
A damped cosine, measured with noise that grows over time — so the input variances are not all equal. That is the interesting case: a single “average” error bar would not capture it.
The input is a one-dimensional scipp.DataArray whose data has variances and whose evenly spaced coordinate carries a unit. A plain DataArray’s variances are treated as uncorrelated; here we build a CovDataArray instead, which carries a full covariance matrix (diagonal in this case, but it need not be).
[2]:
rng = np.random.default_rng(42)
N, dt = 128, 0.05
t = np.arange(N) * dt
truth = np.exp(-t / 1.5) * np.cos(2 * np.pi * 1.5 * t)
sigma = 0.05 + 0.15 * (t / t.max()) # noise grows with time
measured = truth + rng.normal(0.0, sigma)
[3]:
data = CovDataArray(
data=CovVariable(dims=["t"], values=measured, covariance=np.diag(sigma**2)),
coords={"t": sc.array(dims=["t"], values=t, unit="s")},
)
data
[3]:
- t: 128
- t(t)float64s0.0, 0.05, ..., 6.300, 6.350
Values:
array([0. , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 , 0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1. , 1.05, 1.1 , 1.15, 1.2 , 1.25, 1.3 , 1.35, 1.4 , 1.45, 1.5 , 1.55, 1.6 , 1.65, 1.7 , 1.75, 1.8 , 1.85, 1.9 , 1.95, 2. , 2.05, 2.1 , 2.15, 2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 , 2.55, 2.6 , 2.65, 2.7 , 2.75, 2.8 , 2.85, 2.9 , 2.95, 3. , 3.05, 3.1 , 3.15, 3.2 , 3.25, 3.3 , 3.35, 3.4 , 3.45, 3.5 , 3.55, 3.6 , 3.65, 3.7 , 3.75, 3.8 , 3.85, 3.9 , 3.95, 4. , 4.05, 4.1 , 4.15, 4.2 , 4.25, 4.3 , 4.35, 4.4 , 4.45, 4.5 , 4.55, 4.6 , 4.65, 4.7 , 4.75, 4.8 , 4.85, 4.9 , 4.95, 5. , 5.05, 5.1 , 5.15, 5.2 , 5.25, 5.3 , 5.35, 5.4 , 5.45, 5.5 , 5.55, 5.6 , 5.65, 5.7 , 5.75, 5.8 , 5.85, 5.9 , 5.95, 6. , 6.05, 6.1 , 6.15, 6.2 , 6.25, 6.3 , 6.35])
- (t)float64𝟙1.015, 0.809, ..., 0.277, -0.252σ = 0.05, 0.051, ..., 0.199, 0.2
Values:
array([ 1.01523585, 0.80856832, 0.58917249, 0.1919087 , -0.37721268, -0.67135198, -0.77136124, -0.80056683, -0.62064785, -0.38804442, 0.05435649, 0.36362863, 0.54653771, 0.71403224, 0.6275031 , 0.37069363, 0.20668971, -0.15596034, -0.25998516, -0.47658027, -0.52702706, -0.49339642, -0.18942248, -0.08459736, 0.10529232, 0.27930302, 0.44273893, 0.43149024, 0.35242448, 0.20897224, 0.18296753, -0.19673971, -0.32339875, -0.40117944, -0.25066536, -0.1170766 , -0.10361651, -0.03315085, 0.0873915 , 0.3053255 , 0.33587422, 0.28062699, 0.0786571 , 0.06071052, -0.05938984, -0.13521921, -0.11433556, -0.18257792, -0.09090245, -0.08136413, 0.03152995, 0.15252741, -0.01940359, 0.13280315, 0.10368981, 0.03960625, 0.01583035, 0.15199307, -0.18763645, -0.00878342, -0.33873726, -0.15750318, -0.054362 , 0.05377519, 0.12592388, 0.18157928, 0.06075973, 0.04614588, 0.19566747, 0.0203608 , -0.16925445, -0.19428093, -0.19755462, -0.01894012, -0.06114406, 0.03764645, -0.08424852, 0.03435813, 0.13256956, 0.01967711, 0.13548208, -0.03654173, -0.01510723, -0.046674 , -0.19722568, 0.03164807, -0.1252535 , -0.05243724, 0.03094935, 0.04589609, 0.10399917, 0.00635232, -0.02948021, 0.03175221, -0.23026097, -0.20492828, -0.20351416, -0.17028141, 0.04384758, -0.18401398, -0.09924699, 0.1892044 , -0.0803495 , 0.12154759, -0.15171368, -0.0143966 , -0.13866299, -0.0318959 , 0.17130967, -0.29674349, 0.07816205, 0.03183023, -0.12764973, -0.28814465, -0.00795747, -0.11369385, 0.03704517, 0.00727887, 0.31483683, -0.028737 , -0.17792187, 0.05036901, 0.0527715 , 0.26800869, 0.15910985, 0.05956833, 0.27667063, -0.25207801])
Variances (σ²):
array([0.0025 , 0.00261951, 0.0027418 , 0.00286689, 0.00299476, 0.00312543, 0.00325888, 0.00339513, 0.00353416, 0.00367599, 0.0038206 , 0.00396801, 0.0041182 , 0.00427119, 0.00442696, 0.00458553, 0.00474688, 0.00491103, 0.00507797, 0.00524769, 0.00542021, 0.00559551, 0.00577361, 0.00595449, 0.00613817, 0.00632463, 0.00651389, 0.00670593, 0.00690077, 0.00709839, 0.00729881, 0.00750202, 0.00770801, 0.0079168 , 0.00812837, 0.00834274, 0.00855989, 0.00877984, 0.00900257, 0.0092281 , 0.00945641, 0.00968752, 0.00992141, 0.0101581 , 0.01039758, 0.01063984, 0.0108849 , 0.01113274, 0.01138338, 0.0116368 , 0.01189302, 0.01215202, 0.01241382, 0.01267841, 0.01294578, 0.01321595, 0.0134889 , 0.01376465, 0.01404318, 0.01432451, 0.01460862, 0.01489553, 0.01518523, 0.01547771, 0.01577299, 0.01607105, 0.01637191, 0.01667555, 0.01698199, 0.01729121, 0.01760323, 0.01791804, 0.01823563, 0.01855602, 0.01887919, 0.01920516, 0.01953391, 0.01986546, 0.0201998 , 0.02053692, 0.02087684, 0.02121954, 0.02156504, 0.02191332, 0.0222644 , 0.02261827, 0.02297492, 0.02333437, 0.0236966 , 0.02406163, 0.02442944, 0.02480005, 0.02517345, 0.02554963, 0.02592861, 0.02631037, 0.02669493, 0.02708227, 0.02747241, 0.02786534, 0.02826105, 0.02865956, 0.02906085, 0.02946494, 0.02987181, 0.03028148, 0.03069394, 0.03110918, 0.03152722, 0.03194804, 0.03237166, 0.03279807, 0.03322726, 0.03365925, 0.03409402, 0.03453159, 0.03497194, 0.03541509, 0.03586103, 0.03630975, 0.03676127, 0.03721557, 0.03767267, 0.03813256, 0.03859523, 0.0390607 , 0.03952895, 0.04 ])
[4]:
fig, ax = plt.subplots(figsize=(8, 3))
ax.errorbar(t, measured, yerr=sigma, fmt="o", ms=3, lw=1, alpha=0.7, label="measured")
ax.plot(t, truth, "k-", lw=1.5, label="truth")
ax.set(xlabel="t / s", ylabel="signal")
ax.legend()
plt.show()
Exact propagation#
The DFT is linear, so an input covariance maps onto the output exactly: cos @ cov @ cos.T for the real part and sin @ cov @ sin.T for the imaginary part. No sampling, no tuning parameter — just one matrix product.
The result is a scipp.DataGroup with real and imag entries on an omega axis of angular frequencies (here in rad/s, from the s coordinate). The values are (2 / N) * np.fft.fft(values), so a cosine of amplitude A shows up with height A. The transform correlates different frequencies, so each entry is a CovDataArray with the full covariance matrix in .covariance; .variances is its diagonal.
[5]:
ana = analytical.dft(data, "t")
[6]:
def plot_spectrum(ds, ax, label, **kw):
"""Plot the positive-frequency real part with a 1-sigma band."""
w = ds.coords["omega"].values
keep = w >= 0
v = ds.values[keep]
e = np.sqrt(ds.variances[keep])
ax.plot(w[keep], v, label=label, **kw)
ax.fill_between(w[keep], v - e, v + e, alpha=0.25, lw=0, **kw)
fig, ax = plt.subplots(figsize=(8, 3))
plot_spectrum(ana["real"], ax, "real", color="C0")
plot_spectrum(ana["imag"], ax, "imag", color="C1")
ax.set(xlabel="omega / rad s$^{-1}$", ylabel="F")
ax.legend()
plt.show()
Monte Carlo#
The numerical route draws n_samples realisations from a multivariate normal built from the input values and covariance, transforms every one of them, and takes the mean and covariance of the results.
It is (usually) slower and stochastic, but it makes no assumptions about how the uncertainty propagates — which is what makes it the reference when you want to check the analytical result, or when the analytical assumptions do not hold.
[7]:
num = numerical.dft(data, "t", n_samples=5_000)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))
plot_spectrum(ana["real"], axes[0], "analytical", color="C0")
plot_spectrum(num["real"], axes[0], "numerical", color="C1")
axes[0].set(xlabel="omega / rad s$^{-1}$", ylabel="Re[F]")
axes[0].legend()
w = ana["real"].coords["omega"].values
keep = w >= 0
axes[1].plot(
w[keep], np.sqrt(ana["real"].variances[keep]), "o-", ms=3, label="analytical"
)
axes[1].plot(
w[keep], np.sqrt(num["real"].variances[keep]), "x-", ms=4, label="numerical"
)
axes[1].set(xlabel="omega / rad s$^{-1}$", ylabel="sigma(Re[F])")
axes[1].legend()
plt.tight_layout()
plt.show()
[8]:
sd_a = np.sqrt(ana["real"].variances)
sd_n = np.sqrt(num["real"].variances)
print(f"max relative difference in sigma: {np.abs(sd_n / sd_a - 1).max():.2%}")
max relative difference in sigma: 3.01%
The two agree to within Monte Carlo noise, and the sampling error shrinks as n_samples grows — the usual \(1/\sqrt{n}\) convergence.
[9]:
for n in (100, 500, 2_500, 12_500):
sd = np.sqrt(numerical.dft(data, "t", n_samples=n)["real"].variances)
print(
f"n_samples = {n:6d} max relative difference = {np.abs(sd / sd_a - 1).max():.2%}"
)
n_samples = 100 max relative difference = 14.90%
n_samples = 500 max relative difference = 9.76%
n_samples = 2500 max relative difference = 2.82%
n_samples = 12500 max relative difference = 1.56%
Fast Fourier transform#
numerical.fft samples in exactly the same way but uses np.fft.fft instead of an explicit DFT matrix, which is the difference between O(N²) and O(N log N) per sample.
[10]:
fast = numerical.fft(data, "t", n_samples=5_000)
fig, ax = plt.subplots(figsize=(8, 3))
plot_spectrum(fast["real"], ax, "fft", color="C2")
ax.set(xlabel="omega / rad s$^{-1}$", ylabel="Re[F]")
ax.legend()
plt.show()
Likewise, analytical.fft gives exactly the same result as analytical.dft, but computes the propagated covariance with FFTs, which scales as O(N² log N) rather than with dense matrix products.