NFFT的傅立葉係數 - 非均勻快速傅里葉變換？

Question

我想在python 2.7中使用包pynfft來做非均勻快速傅里葉變換（nfft）。我只學了兩個月python，所以我有一些困難。

這是我的代碼：

import numpy as np 
from pynfft.nfft import NFFT 

#loading data, 104 lines 
t_diff, x_diff = np.loadtxt('data/analysis/amplitudes.dat', unpack = True) 

N = [13,8] 
M = 52 

#fourier coefficients 
f_hat = np.fft.fft(x_diff)/(2*M) 

#instantiation 
plan = NFFT(N,M) 

#precomputation 
x = t_diff 
plan.x = x 
plan.precompute() 

# vector of non uniform samples 
f = x_diff[0:M] 

#execution 
plan.f = f 
plan.f_hat = f_hat 
f = plan.trafo() 

我基本上跟隨我在pynfft教程（http://pythonhosted.org/pyNFFT/tutorial.html）中的說明。

我需要nfft，因爲我的數據採取的時間間隔不是常數（我的意思是，第一個測量是在t，第二個在dt之後，第三個在dt + dt'之後，dt'不同於dt等）。

pynfft包在執行之前需要傅立葉係數（「f_hat」）的向量，所以我使用numpy.fft來計算它，但我不確定這個過程是否正確。有沒有另一種方式來做到這一點（也許與nfft）？

我也想計算頻率;我知道用numpy.fft有一個命令：對於pynfft是否也有類似的情況？我在教程中沒有找到任何東西。

謝謝你的任何建議，你可以給我。

Answer 1

下面是一個工作示例，從here採取：

首先我們定義我們要重建的功能，這是四個諧波的總和：我們初始化

import numpy as np 
import matplotlib.pyplot as plt 

np.random.seed(12345) 

%pylab inline --no-import-all 

# function we want to reconstruct 
k=[1,5,10,30] # modulating coefficients 
def myf(x,k): 
    return sum(np.sin(x*k0*(2*np.pi)) for k0 in k) 

x=np.linspace(-0.5,0.5,1000) # 'continuous' time/spatial domain; -0.5<x<+0.5 
y=myf(x,k)      # 'true' underlying trigonometric function 

fig=plt.figure(1,(20,5)) 
ax =fig.add_subplot(111) 

ax.plot(x,y,'red') 
ax.plot(x,y,'r.') 

         # we should sample at a rate of >2*~max(k) 
M=256     # number of nodes 
N=128     # number of Fourier coefficients 

nodes =np.random.rand(M)-0.5 # non-uniform oversampling 
values=myf(nodes,k)  # nodes&values will be used below to reconstruct 
         # original function using the Solver 

ax.plot(nodes,values,'bo') 

ax.set_xlim(-0.5,+0.5) 

的和運行規劃求解：

from pynfft import NFFT, Solver 

f  = np.empty(M,  dtype=np.complex128) 
f_hat = np.empty([N,N], dtype=np.complex128) 

this_nfft = NFFT(N=[N,N], M=M) 
this_nfft.x = np.array([[node_i,0.] for node_i in nodes]) 
this_nfft.precompute() 

this_nfft.f = f 
ret2=this_nfft.adjoint() 

print this_nfft.M # number of nodes, complex typed 
print this_nfft.N # number of Fourier coefficients, complex typed 
#print this_nfft.x # nodes in [-0.5, 0.5), float typed 


this_solver = Solver(this_nfft) 
this_solver.y = values   # '''right hand side, samples.''' 

#this_solver.f_hat_iter = f_hat # assign arbitrary initial solution guess, default is 0 

this_solver.before_loop()  # initialize solver internals 

while not np.all(this_solver.r_iter < 1e-2): 
this_solver.loop_one_step() 

最後，我們顯示頻率：

import matplotlib.pyplot as plt 

fig=plt.figure(1,(20,5)) 
ax =fig.add_subplot(111) 


foo=[ np.abs(this_solver.f_hat_iter[i][0])**2 for i in range(len(this_solver.f_hat_iter)) ] 

ax.plot(np.abs(np.arange(-N/2,+N/2,1)),foo) 

歡呼聲