Calcolare l'incertezza in FFT ampiezza

Question

mio Python problema di programmazione è la seguente:

Voglio creare una serie di risultati di misura. Ogni risultato può essere descritto come una distribuzione normale per cui il valore medio è il risultato della misurazione stessa e la deviazione standard è la sua incertezza.

pseudo codice potrebbe essere:

x1 = N(result1, unc1) 
x2 = N(result2, unc2) 
... 

x = array(x1, x2, ..., xN) 

di quanto vorrei per calcolare la FFT di x:

f = numpy.fft.fft(x) 

Quello che voglio è che l'incertezza delle misure contenute in x si propaga attraverso il calcolo FFT in modo che f sia una matrice di ampiezze con la loro incertezza come questa:

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

Potete suggerirmi un modo per farlo?

Answer 1

Ogni coefficiente di Fourier calcolato dalla trasformata discreta di Fourier dell'array x è una combinazione lineare degli elementi di x; vedi la formula per X_k sul wikipedia page on the discrete Fourier transform, che scriverò come

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

(Cioè, X è la trasformata di Fourier discreta x.) Se x_n distribuzione normale con mu_n media e varianza sigma_n ** 2, quindi un po 'di algebra mostra che la varianza X_k è la somma delle varianze di x_n

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

In altre parole, la varianza è uguale per ogni Fouri er coefficent; è la somma delle varianze delle misure in x.

Utilizzando il notazione, dove unc(z) è la deviazione standard di z,

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

(Si noti che la distribuzione del grandezza del X_k è il Rice distribution.)

Ecco uno script che dimostra questo risultato In questo esempio, la deviazione standard dei valori x aumenta linearmente da 0,01 a 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() 

La stampa è

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] 

Come previsto, le varianze campionarie dei coefficienti di Fourier sono tutti (circa) uguale alla somma delle varianze misurazione.

Ecco la trama generata dalla sceneggiatura. I diamanti neri sono i coefficienti di Fourier di un singolo vettore x.I punti blu sono i coefficienti di Fourier di 800 realizzazioni di x + noise. È possibile vedere che le nuvole di punti attorno a ciascun coefficiente di Fourier sono approssimativamente simmetriche e tutte le stesse "dimensioni" (eccetto, ovviamente, per i veri coefficienti, che appaiono in questo grafico come linee orizzontali sull'asse reale).