Berechnen Unsicherheit in FFT Amplitude

Question

Mein Python Programmierung Problem ist folgendes:

Ich möchte eine Reihe von Messergebnissen erstellen. Jedes Ergebnis kann als Normalverteilung beschrieben werden, für die der Mittelwert das Messergebnis selbst und die Standardabweichung seine Unsicherheit ist.

Pseudo-Code könnte sein:

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

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

als ich die FFT von x berechnen möchte:

f = numpy.fft.fft(x) 

Was ich will, ist, dass die Unsicherheit der in x enthaltenen Messungen propagiert durch die FFT-Berechnung, so dass f ein Array von Amplituden mit ihrer Unsicherheit ist:

Können Sie mir einen Weg vorschlagen, dies zu tun?

Answer 1

Jeder Fourier-Koeffizient, der durch die diskrete Fourier-Transformation des Arrays x berechnet wird, ist eine Linearkombination der Elemente x; siehe die Formel für x_K auf der wikipedia page on the discrete Fourier transform, die ich werde schreiben als

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

(Das heißt, X die diskrete Fourier-Transformation von x.) Wenn x_n normalerweise mit einem Mittelwert und einer Varianz mu_n sigma_n verteilt wird, ** 2, dann ein wenig Algebra zeigt, dass die Varianz von x_k die Summe der Varianzen von x_n ist

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

Mit anderen Worten, ist die Varianz das gleiche für jede Fouri er Coefficent; es ist die Summe der Varianzen der Messungen in x.

Ihre Schreibweise verwenden, wo unc(z) die Standardabweichung der z ist,

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

(Beachten Sie, dass die Verteilung der Größenordnung von x_k die Rice distribution ist.)

Hier ist ein Skript, das zeigt, dieses Ergebnis. In diesem Beispiel steigt der Standard Abweichung der x Werte linear von 0,01 bis 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() 

der gedruckten Ausgabe ist

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] 

Wie erwartet, werden die Proben Varianzen der Fourier-Koeffizienten sind alle (etwa) gleich der Summe der Messabweichungen.

Hier ist das Diagramm, das vom Skript generiert wird. Die schwarzen Diamanten sind die Fourier-Koeffizienten eines einzelnen Vektors x.Die blauen Punkte sind die Fourier-Koeffizienten von 800 Realisierungen von x + noise. Sie können sehen, dass die Punktwolken um jeden Fourier-Koeffizienten ungefähr symmetrisch sind und alle die gleiche "Größe" (außer natürlich für die echten Koeffizienten, , die in diesem Diagramm als horizontale Linien auf der realen Achse angezeigt werden).