Usage#
fourierror propagates measurement uncertainty through a discrete Fourier transform. It offers three 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 |
one matrix product |
|
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
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 scipp.DataArray whose data has variances and whose coordinate carries a unit.
[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)
data = sc.DataArray(
data=sc.array(dims=["t"], values=measured, variances=sigma**2),
coords={"t": sc.array(dims=["t"], values=t, unit="s")},
)
data
[2]:
scipp.DataArray (3.86 KB)
- 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 ])
[3]:
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.Dataset with real and imag entries, each carrying its own variances, on an omega axis.
[4]:
ana = analytical.dft(data, "t")
ana
[4]:
scipp.Dataset (6.67 KB)
- omega: 128
- omega(omega)float64Hz0.0, 0.982, ..., -1.963, -0.982
Values:
array([ 0. , 0.9817477 , 1.96349541, 2.94524311, 3.92699082, 4.90873852, 5.89048623, 6.87223393, 7.85398163, 8.83572934, 9.81747704, 10.79922475, 11.78097245, 12.76272016, 13.74446786, 14.72621556, 15.70796327, 16.68971097, 17.67145868, 18.65320638, 19.63495408, 20.61670179, 21.59844949, 22.5801972 , 23.5619449 , 24.54369261, 25.52544031, 26.50718801, 27.48893572, 28.47068342, 29.45243113, 30.43417883, 31.41592654, 32.39767424, 33.37942194, 34.36116965, 35.34291735, 36.32466506, 37.30641276, 38.28816047, 39.26990817, 40.25165587, 41.23340358, 42.21515128, 43.19689899, 44.17864669, 45.1603944 , 46.1421421 , 47.1238898 , 48.10563751, 49.08738521, 50.06913292, 51.05088062, 52.03262833, 53.01437603, 53.99612373, 54.97787144, 55.95961914, 56.94136685, 57.92311455, 58.90486225, 59.88660996, 60.86835766, 61.85010537, -62.83185307, -61.85010537, -60.86835766, -59.88660996, -58.90486225, -57.92311455, -56.94136685, -55.95961914, -54.97787144, -53.99612373, -53.01437603, -52.03262833, -51.05088062, -50.06913292, -49.08738521, -48.10563751, -47.1238898 , -46.1421421 , -45.1603944 , -44.17864669, -43.19689899, -42.21515128, -41.23340358, -40.25165587, -39.26990817, -38.28816047, -37.30641276, -36.32466506, -35.34291735, -34.36116965, -33.37942194, -32.39767424, -31.41592654, -30.43417883, -29.45243113, -28.47068342, -27.48893572, -26.50718801, -25.52544031, -24.54369261, -23.5619449 , -22.5801972 , -21.59844949, -20.61670179, -19.63495408, -18.65320638, -17.67145868, -16.68971097, -15.70796327, -14.72621556, -13.74446786, -12.76272016, -11.78097245, -10.79922475, -9.81747704, -8.83572934, -7.85398163, -6.87223393, -5.89048623, -4.90873852, -3.92699082, -2.94524311, -1.96349541, -0.9817477 ])
- imag(omega)float64𝟙0.0, -0.024, ..., -0.007, 0.024σ = 0.0, 0.016, ..., 0.017, 0.016
Values:
array([ 0.00000000e+00, -2.42925147e-02, 7.37462838e-03, 2.77815430e-02, 2.76505479e-02, 3.12804638e-02, 4.76239324e-02, 7.12011227e-02, 9.68103674e-02, 9.22395917e-02, -9.72072332e-02, -7.41970465e-02, -8.28404732e-02, -4.63209916e-02, -4.38481221e-02, -1.22799039e-02, -3.19985584e-02, -2.13160033e-02, -2.96641985e-02, -4.06257881e-02, -3.17651843e-02, -3.22528574e-02, -1.09772966e-02, -5.07539748e-02, -1.26344654e-02, -1.37505555e-02, -3.73460073e-03, -1.58720648e-02, -6.45967210e-03, -1.01083472e-02, -1.00863624e-02, -4.76137106e-04, -5.63879972e-03, 1.59355169e-03, -1.04705132e-02, -1.59965004e-02, -3.29952277e-02, -1.16480707e-02, -5.49759290e-03, 6.63451509e-03, 2.29316465e-03, -9.15515531e-03, -2.07503792e-02, -9.89378473e-03, -1.33797594e-02, 2.22176624e-03, -1.25902096e-02, -5.55705149e-04, -1.73897171e-02, -1.48574012e-02, -4.07827332e-02, -1.74806208e-02, 6.83935236e-03, -1.40633381e-02, -1.20596750e-03, -1.53921595e-02, -1.29611289e-03, 1.72316958e-03, 1.93167602e-02, -4.61407396e-03, -3.61461822e-03, -2.09518991e-02, 8.44340597e-04, 2.91968896e-03, -1.77308425e-16, -2.91968896e-03, -8.44340597e-04, 2.09518991e-02, 3.61461822e-03, 4.61407396e-03, -1.93167602e-02, -1.72316958e-03, 1.29611289e-03, 1.53921595e-02, 1.20596750e-03, 1.40633381e-02, -6.83935236e-03, 1.74806208e-02, 4.07827332e-02, 1.48574012e-02, 1.73897171e-02, 5.55705149e-04, 1.25902096e-02, -2.22176624e-03, 1.33797594e-02, 9.89378473e-03, 2.07503792e-02, 9.15515531e-03, -2.29316465e-03, -6.63451509e-03, 5.49759290e-03, 1.16480707e-02, 3.29952277e-02, 1.59965004e-02, 1.04705132e-02, -1.59355169e-03, 5.63879972e-03, 4.76137106e-04, 1.00863624e-02, 1.01083472e-02, 6.45967210e-03, 1.58720648e-02, 3.73460073e-03, 1.37505555e-02, 1.26344654e-02, 5.07539748e-02, 1.09772966e-02, 3.22528574e-02, 3.17651843e-02, 4.06257881e-02, 2.96641985e-02, 2.13160033e-02, 3.19985584e-02, 1.22799039e-02, 4.38481221e-02, 4.63209916e-02, 8.28404732e-02, 7.41970465e-02, 9.72072332e-02, -9.22395917e-02, -9.68103674e-02, -7.12011227e-02, -4.76239324e-02, -3.12804638e-02, -2.76505479e-02, -2.77815430e-02, -7.37462838e-03, 2.42925147e-02])
Variances (σ²):
array([0.00000000e+00, 2.71689986e-04, 2.75082221e-04, 2.75710404e-04, 2.75930259e-04, 2.76032011e-04, 2.76087272e-04, 2.76120581e-04, 2.76142188e-04, 2.76156989e-04, 2.76167562e-04, 2.76175372e-04, 2.76181298e-04, 2.76185895e-04, 2.76189527e-04, 2.76192442e-04, 2.76194810e-04, 2.76196756e-04, 2.76198368e-04, 2.76199714e-04, 2.76200843e-04, 2.76201793e-04, 2.76202595e-04, 2.76203271e-04, 2.76203839e-04, 2.76204313e-04, 2.76204706e-04, 2.76205025e-04, 2.76205277e-04, 2.76205469e-04, 2.76205603e-04, 2.76205682e-04, 2.76205709e-04, 2.76205682e-04, 2.76205603e-04, 2.76205469e-04, 2.76205277e-04, 2.76205025e-04, 2.76204706e-04, 2.76204313e-04, 2.76203839e-04, 2.76203271e-04, 2.76202595e-04, 2.76201793e-04, 2.76200843e-04, 2.76199714e-04, 2.76198368e-04, 2.76196756e-04, 2.76194810e-04, 2.76192442e-04, 2.76189527e-04, 2.76185895e-04, 2.76181298e-04, 2.76175372e-04, 2.76167562e-04, 2.76156989e-04, 2.76142188e-04, 2.76120581e-04, 2.76087272e-04, 2.76032011e-04, 2.75930259e-04, 2.75710404e-04, 2.75082221e-04, 2.71689986e-04, 1.81171480e-31, 2.71689986e-04, 2.75082221e-04, 2.75710404e-04, 2.75930259e-04, 2.76032011e-04, 2.76087272e-04, 2.76120581e-04, 2.76142188e-04, 2.76156989e-04, 2.76167562e-04, 2.76175372e-04, 2.76181298e-04, 2.76185895e-04, 2.76189527e-04, 2.76192442e-04, 2.76194810e-04, 2.76196756e-04, 2.76198368e-04, 2.76199714e-04, 2.76200843e-04, 2.76201793e-04, 2.76202595e-04, 2.76203271e-04, 2.76203839e-04, 2.76204313e-04, 2.76204706e-04, 2.76205025e-04, 2.76205277e-04, 2.76205469e-04, 2.76205603e-04, 2.76205682e-04, 2.76205709e-04, 2.76205682e-04, 2.76205603e-04, 2.76205469e-04, 2.76205277e-04, 2.76205025e-04, 2.76204706e-04, 2.76204313e-04, 2.76203839e-04, 2.76203271e-04, 2.76202595e-04, 2.76201793e-04, 2.76200843e-04, 2.76199714e-04, 2.76198368e-04, 2.76196756e-04, 2.76194810e-04, 2.76192442e-04, 2.76189527e-04, 2.76185895e-04, 2.76181298e-04, 2.76175372e-04, 2.76167562e-04, 2.76156989e-04, 2.76142188e-04, 2.76120581e-04, 2.76087272e-04, 2.76032011e-04, 2.75930259e-04, 2.75710404e-04, 2.75082221e-04, 2.71689986e-04]) - real(omega)float64𝟙-0.004, 0.017, ..., 0.026, 0.017σ = 0.023, 0.017, ..., 0.017, 0.017
Values:
array([-4.34913327e-03, 1.70044689e-02, 2.61698261e-02, 2.47495808e-02, 1.21130717e-02, 1.58014536e-02, 2.74687203e-02, 3.75093372e-03, 3.40242829e-02, 1.35808242e-01, 2.11891923e-01, 3.67639399e-02, 1.54801213e-02, 2.87611653e-02, 2.39673096e-02, -8.63325519e-03, -4.36247226e-03, -1.60338161e-02, -2.51248172e-04, -1.47334190e-03, 5.85973102e-04, 6.69438609e-03, 1.80919304e-02, 8.68461768e-03, 9.01966344e-03, 1.06860341e-02, 2.10825463e-02, -1.55139779e-02, 3.72799128e-03, 1.41298161e-02, 1.54831104e-02, -1.31056705e-02, -1.40841388e-02, 6.48144273e-04, -5.44669579e-03, -4.82316152e-03, 1.16764933e-02, 1.33493433e-02, 2.64095015e-02, 1.90241199e-02, -1.45994720e-04, 3.02035304e-03, -2.66816509e-03, 1.84130850e-02, -7.27181377e-04, 1.25843150e-02, 2.44027451e-03, 1.43517773e-02, 1.84708667e-02, -1.04645671e-02, 2.22793188e-02, 4.29117503e-02, 2.01287398e-02, 8.94655554e-03, 9.00253056e-03, 8.69989906e-03, 4.44619373e-02, 1.55398008e-02, 1.21964807e-02, 5.44681874e-03, -1.65680772e-03, 9.49447943e-03, 8.90079732e-03, 3.94386994e-02, 1.40474515e-02, 3.94386994e-02, 8.90079732e-03, 9.49447943e-03, -1.65680772e-03, 5.44681874e-03, 1.21964807e-02, 1.55398008e-02, 4.44619373e-02, 8.69989906e-03, 9.00253056e-03, 8.94655554e-03, 2.01287398e-02, 4.29117503e-02, 2.22793188e-02, -1.04645671e-02, 1.84708667e-02, 1.43517773e-02, 2.44027451e-03, 1.25843150e-02, -7.27181377e-04, 1.84130850e-02, -2.66816509e-03, 3.02035304e-03, -1.45994720e-04, 1.90241199e-02, 2.64095015e-02, 1.33493433e-02, 1.16764933e-02, -4.82316152e-03, -5.44669579e-03, 6.48144273e-04, -1.40841388e-02, -1.31056705e-02, 1.54831104e-02, 1.41298161e-02, 3.72799128e-03, -1.55139779e-02, 2.10825463e-02, 1.06860341e-02, 9.01966344e-03, 8.68461768e-03, 1.80919304e-02, 6.69438609e-03, 5.85973102e-04, -1.47334190e-03, -2.51248172e-04, -1.60338161e-02, -4.36247226e-03, -8.63325519e-03, 2.39673096e-02, 2.87611653e-02, 1.54801213e-02, 3.67639399e-02, 2.11891923e-01, 1.35808242e-01, 3.40242829e-02, 3.75093372e-03, 2.74687203e-02, 1.58014536e-02, 1.21130717e-02, 2.47495808e-02, 2.61698261e-02, 1.70044689e-02])
Variances (σ²):
array([0.0005478 , 0.00027611, 0.00027272, 0.00027209, 0.00027187, 0.00027177, 0.00027171, 0.00027168, 0.00027166, 0.00027164, 0.00027163, 0.00027162, 0.00027162, 0.00027161, 0.00027161, 0.00027161, 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.00027161, 0.00027161, 0.00027161, 0.00027162, 0.00027162, 0.00027163, 0.00027164, 0.00027166, 0.00027168, 0.00027171, 0.00027177, 0.00027187, 0.00027209, 0.00027272, 0.00027611, 0.0005478 , 0.00027611, 0.00027272, 0.00027209, 0.00027187, 0.00027177, 0.00027171, 0.00027168, 0.00027166, 0.00027164, 0.00027163, 0.00027162, 0.00027162, 0.00027161, 0.00027161, 0.00027161, 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.00027159, 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.0002716 , 0.00027161, 0.00027161, 0.00027161, 0.00027162, 0.00027162, 0.00027163, 0.00027164, 0.00027166, 0.00027168, 0.00027171, 0.00027177, 0.00027187, 0.00027209, 0.00027272, 0.00027611])
[5]:
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["real"].values[keep]
e = np.sqrt(ds["real"].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, ax, "analytical", color="C0")
ax.set(xlabel="omega / rad s$^{-1}$", ylabel="Re[F]")
ax.legend()
plt.show()
Monte Carlo#
The numerical route draws n_samples realisations from a multivariate normal built from the input values and variances, 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.
[10]:
num = numerical.dft(data, "t", n_samples=5_000)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))
plot_spectrum(ana, axes[0], "analytical", color="C0")
plot_spectrum(num, axes[0], "numerical", color="C1")
axes[0].set(xlabel="omega / rad s$^{-1}$", ylabel="Re[F]")
axes[0].legend()
w = ana.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()
[11]:
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: 2.44%
The two agree to within Monte Carlo noise, and the sampling error shrinks as n_samples grows — the usual \(1/\sqrt{n}\) convergence.
[12]:
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 = 15.38%
n_samples = 500 max relative difference = 9.37%
n_samples = 2500 max relative difference = 3.45%
n_samples = 12500 max relative difference = 1.58%
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.
We are in the process of developing an FFT analytical form.
[15]:
fast = numerical.fft(data, "t", n_samples=5_000)
fig, ax = plt.subplots(figsize=(8, 3))
plot_spectrum(fast, ax, "fft", color="C2")
ax.set(xlabel="omega / rad s$^{-1}$", ylabel="Re[F]")
ax.legend()
plt.show()
