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的命令？我在教程中没有找到任何信息。

非常感谢您能给我的任何建议。

- elle

Answer 1

这里是一个工作示例，取自这里：

首先，我们定义我们想要重构的函数，它是四个谐波的和：

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)

干杯

- Vincenzo Lavorini