My Python programming problem is the following:

I want to create an array of measurement results. Each result can be described as a normal distribution for which the mean value is the measurement result itself and the standard deviation is its uncertainty.

Pseudo code could be:

```
x1 = N(result1, unc1)
x2 = N(result2, unc2)
...
x = array(x1, x2, ..., xN)
```

Than I would like to calculate the FFT of x:

```
f = numpy.fft.fft(x)
```

What I want is that the uncertainty of the measurements contained in x is propagated through the FFT calculation so that f is an array of amplitudes along with their uncertainty like this:

```
f = (a +/- unc(a), b +/- unc(b), ...)
```

Can you suggest me a way to do this?

pythonnumpyfftuncertainty answered 4 years ago Warren Weckesser #1

Each Fourier coefficient computed by the discrete Fourier transform
of the array `x`

is a linear combination of the elements of `x`

; see
the formula for X_k on the wikipedia page on the discrete Fourier transform,
which I'll write as

```
X_k = sum_(n=0)^(n=N-1) [ x_n * exp(-i*2*pi*k*n/N) ]
```

(That is, `X`

is the discrete Fourier transform of `x`

.)
If x_n is normally distributed with mean mu_n and variance sigma_n**2,
then a little bit of algebra shows that the variance of X_k is the sum
of the variances of x_n

```
Var(X_k) = sum_(n=0)^(n=N-1) sigma_n**2
```

In other words, the variance is the same for each Fourier coefficent;
it is the sum of the variances of the measurements in `x`

.

Using your notation, where `unc(z)`

is the standard deviation of `z`

,

```
unc(X_0) = unc(X_1) = ... = unc(X_(N-1)) = sqrt(unc(x1)**2 + unc(x2)**2 + ...)
```

(Note that the distribution of the *magnitude* of X_k is the Rice distribution.)

Here's a script that demonstrates this result. In this example, the standard
deviation of the `x`

values increase linearly from 0.01 to 0.5.

```
import numpy as np
from numpy.fft import fft
import matplotlib.pyplot as plt
np.random.seed(12345)
n = 16
# Create 'x', the vector of measured values.
t = np.linspace(0, 1, n)
x = 0.25*t - 0.2*t**2 + 1.25*np.cos(3*np.pi*t) + 0.8*np.cos(7*np.pi*t)
x[:n//3] += 3.0
x[::4] -= 0.25
x[::3] += 0.2
# Compute the Fourier transform of x.
f = fft(x)
num_samples = 5000000
# Suppose the std. dev. of the 'x' measurements increases linearly
# from 0.01 to 0.5:
sigma = np.linspace(0.01, 0.5, n)
# Generate 'num_samples' arrays of the form 'x + noise', where the standard
# deviation of the noise for each coefficient in 'x' is given by 'sigma'.
xn = x + sigma*np.random.randn(num_samples, n)
fn = fft(xn, axis=-1)
print "Sum of input variances: %8.5f" % (sigma**2).sum()
print
print "Variances of Fourier coefficients:"
np.set_printoptions(precision=5)
print fn.var(axis=0)
# Plot the Fourier coefficient of the first 800 arrays.
num_plot = min(num_samples, 800)
fnf = fn[:num_plot].ravel()
clr = "#4080FF"
plt.plot(fnf.real, fnf.imag, 'o', color=clr, mec=clr, ms=1, alpha=0.3)
plt.plot(f.real, f.imag, 'kD', ms=4)
plt.grid(True)
plt.axis('equal')
plt.title("Fourier Coefficients")
plt.xlabel("$\Re(X_k)$")
plt.ylabel("$\Im(X_k)$")
plt.show()
```

The printed output is

```
Sum of input variances: 1.40322
Variances of Fourier coefficients:
[ 1.40357 1.40288 1.40331 1.40206 1.40231 1.40302 1.40282 1.40358
1.40376 1.40358 1.40282 1.40302 1.40231 1.40206 1.40331 1.40288]
```

As expected, the sample variances of the Fourier coefficients are all (approximately) the same as the sum of the measurement variances.

Here's the plot generated by the script. The black diamonds are the
Fourier coefficients of a single `x`

vector. The blue dots are the
Fourier coefficients of 800 realizations of `x + noise`

. You can see that
the point clouds around each Fourier coefficent are roughly symmetric
and all the same "size" (except, of course, for the real coeffcients,
which show up in this plot as horizontal lines on the real axis).