FFT 진폭의 불확실성을 계산하십시오.

Question

프로그래밍 문제는 다음과 같습니다.

측정 결과 배열을 만들고 싶습니다. 각 결과는 평균값이 측정 결과 자체이고 표준 편차가 불확도 인 정규 분포로 설명 될 수 있습니다.

의사 코드가 될 수있다 :

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

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

내가 x의 FFT을 계산하고 싶습니다보다 :

f = numpy.fft.fft(x) 

내가 원하는 X에 포함 된 측정의 불확실성을 통해 전파된다는 점이다 FFT 계산은 f가 다음과 같이 불확도와 함께 진폭의 배열이되도록합니다.

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

나에게이 방법을 제안 해 줄 수 있습니까?

Answer 1

배열 x의 이산 푸리에 변환 에 의해 계산 된 각 푸리에 계수는 x의 요소의 선형 조합입니다. 제가

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

같이 쓸 겁니다 wikipedia page on the discrete Fourier transform, 에 X_k의 수식 표시 (즉, X 이산 푸리에 x의 변환이다.) x_n 정상적으로 평균 mu_n 및 분산 sigma_n 함께 배포되면 ** 2 후 대수 조금 X_k의 편차가 합계 x_n 말하면

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

을의 차이의 임을 나타내고, 편차 각 Fouri 대해 동일 어 계수기; x에 측정 값의 분산의 합계입니다. unc(z)이 z의 표준 편차 당신의 표기법을, 사용

는

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

합니다 ( 크기의 분포가 X_k의는 Rice distribution 있음을 유의하십시오.) 여기

을 보여줍니다 스크립트입니다 이 결과. 이 예에서는 표준 x 값의 편차가 0.01에서 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() 

인쇄물은 예상 한 바와 같이, 푸리에 계수의 샘플 차이는 측정 차이의 합으로서 전체 (대략) 동일한 이다

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] 

이다.

다음은 스크립트에 의해 생성 된 플롯입니다. 검은 색 다이아몬드는 단일 벡터의 푸리에 계수입니다. x 벡터입니다.파란 점은 의 800 실현의 푸리에 계수입니다. 각 푸리에 coefficent 주위의 점 구름은 대략 대칭입니다 및 모든 "크기"(물론 실제 coeffcients에 대해서는 은 실제 축에 수평선으로 표시됨).