如何解决示例 python nfft 傅立叶变换 - 信号重建归一化问题
我为 nfft 和 scipy.fft 编写了一个完整的工作示例。 在这两种情况下，我都从一个带有一点噪声的简单一维正弦信号开始，进行傅立叶变换，然后倒退并重建原始信号。
这是我的代码，尽我所能：
import numpy
import nfft
import scipy
import scipy.fft
import matplotlib.pyplot as plt
if True: #<--- Ensure non-global namespace
#Define signal:
x = -0.5 + numpy.random.rand(1000)
#x = numpy.linspace(-.5,.5,1000) #--> in case we want to run uniform time domain
f = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( 1000 ) #Add some 'y' randomness to the sample
#prepare wavenumbers for transform:
N = len(x)
k = - N // 2 + numpy.arange(N)
#print ('k',k) #---> Uniform Steps [-500,-499,...0...,499,500]
f_k = nfft.nfft_adjoint(x,f,len(k),truncated=False )
#plot transform
plt.figure()
plt.plot(k,f_k.real,label='real')
plt.plot(k,f_k.imag,label='imag')
plt.legend()
#Reconstruct the original signal with nfft
f_recon = nfft.nfft( x,f_k ) / 2000
#Plot original vs reconstructed
plt.figure()
plt.title('nfft')
plt.scatter(x,label='f(x)')
plt.scatter(x,f_recon,label='f_recon(x)',marker='+')
plt.legend()
if True: #<--- Ensure non-global namespace
#Define signal:
x = numpy.linspace(-.5,1000)
f = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( 1000 ) #Add some 'y' randomness to the sample
#prepare wavenumbers for transform:
N = len(x)
TimeSpacing = x[1] - x[0]
k = scipy.fft.fftfreq(N,TimeSpacing)
#print ('k',k) #---> Confusing steps: [0,1,...500,-500,...-1]
f_k = scipy.fft.fft(f)
#plot transform
plt.figure()
plt.plot(k,label='imag')
plt.legend()
#Reconstruct the original signal with scipy.fft
f_recon = scipy.fft.ifft(f_k)
#Plot original vs reconstructed
plt.figure()
plt.title('scipy.fft')
plt.scatter(x,marker='+')
plt.legend()
plt.show()
以下是相关的生成图：
nfft 重建似乎无法正常化。 我随意将幅度除以 2000 只是为了让它们绘制得很好。 什么是正确的归一化常数？
nfft 似乎也没有重现原点。 即使我得到了正确的归一化常数，我也无法将原始点恢复到这里。
我做错了什么，我该如何解决？
解决方法
上面提到的包没有实现逆nfft
ndft 是 f_hat @ np.exp(-2j * np.pi * x * k[:,None])
ndft_adjoint 是 f @ np.exp(2j * np.pi * k * x[:,None])
让 k = -N//2 + np.arange(N) 和 A = np.exp(-2j * np.pi * k * k[:,None])
A @ np.conj(A) = N * np.eye(N)（以数字方式检查）
因此，对于随机 x，伴随变换等于逆变换。给定的参考文件提供了一些选项，我实现了 Algorithm 1 CGNE,from page 9
import numpy as np # I have the habit to use np
def nfft_inverse(x,y,N,w = 1,L=100):
f = np.zeros(N,dtype=np.complex128);
r = y - nfft.nfft(x,f);
p = nfft.nfft_adjoint(x,r,N);
r_norm = np.sum(abs(r)**2 * w)
for l in range(L):
p_norm = np.sum(abs(p)**2 * w);
alpha = r_norm / p_norm
f += alpha * w * p;
r = y - nfft.nfft(x,f)
r_norm_2 = np.sum(abs(r)**2 * w)
beta = r_norm_2 / r_norm
p = beta * p + nfft.nfft_adjoint(x,w * r,N)
r_norm = r_norm_2;
#print(l,r_norm)
return f;
算法收敛缓慢且效果不佳
plt.figure(figsize=(14,7))
plt.title('inverse nfft error histogram')
#plt.scatter(x,f_hat,label='f(x)')
h_hat = nfft_inverse(x,f,L = 1)
plt.hist(f_hat - numpy.real(h_hat),bins=30,label='1 iteration')
h_hat = nfft_inverse(x,L = 10)
plt.hist(f_hat - numpy.real(h_hat),label='10 iterations')
h_hat = nfft_inverse(x,L = 1000)
plt.hist(f_hat - numpy.real(h_hat),label='1000 iterations')
plt.xlabel('error')
plt.ylabel('occurrencies')
plt.legend()
我也尝试使用 scipy minimization，显式地最小化残差 ||nfft(x,f) - y||**2
import numpy as np # the habit
import scipy.optimize
def nfft_gradient_descent(x,L=10,tol=1e-8,method='CG'):
'''
compute $min || A @ f - y ||**2 via gradient descent
the gradient is
`A^H @ (A @ f - y)`
Multiply by A using nfft.nfft
'''
def cost(fpack):
f = fpack[0::2] + 1j * fpack[1::2]
u = np.sum(np.abs(nfft.nfft(x,f) - y)**2)
return u
def grad(fpack):
f = fpack[0::2] + 1j * fpack[1::2]
r = nfft.nfft(x,f) - y
u = nfft.nfft_adjoint(x,N)
return np.stack([np.real(u),np.imag(u)],axis=1).reshape(-1)
x0 = np.zeros([N,2])
result = scipy.optimize.minimize(cost,x0=x0,jac=grad,tol=tol,method=method,options={'maxiter': L,'disp': True})
return result.x[0::2] + 1j * result.x[1::2];
结果看起来很相似，您可以根据自己的需要尝试不同的方法或参数。但我相信转换是病态的，因为转换后的残差大大减少，但重建值的残差很大。
编辑 1
你发现算法没有真正的逆是基本正确的吗？我无法获得我的原始积分？ x != nfft(nfft_adjoint(x))
请查看参考文献 paper
的第 2.3 节
数值比较
Cris Luengo answer 提到了另一种可能性，即，您可以使用 f 在等距点重建重采样版本，而不是在点 x 处重建 ifft。所以你已经有了三个选项，我会做一个快速的比较。请记住，此处显示的图基于在 16k 个样本中计算的 NFFT，而这里我使用的是 1k 个样本。
由于FFT方法使用的点不同，我们无法与原始信号进行比较，因此我要做的是与没有噪声的谐波函数进行比较。噪声的方差为 0.01，因此精确的重建会导致此均方误差。
N = 1024
x = -0.5 + numpy.random.rand(N)
f_hat = numpy.sin(10 * 2 * numpy.pi * x) + .1*numpy.random.randn( N ) #Add some 'y' randomness to the sample
k = - N // 2 + numpy.arange(N)
f = nfft.nfft(x,f_hat)
print('nfft_inverse')
h_hat = nfft_inverse(x,len(x),L = 10)
print('10 iterations: ',np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_inverse(x,L = 100)
print('100 iterations: ',L = 1000)
print('1000 iterations: ',np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
print('nfft_gradient_descent')
h_hat = nfft_gradient_descent(x,np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x,np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
#Reconstruct the original at a spaced grid based on nfft result using ifft
f_recon = - numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k))) / (N / N2)
x_recon = k / N;
print('using IFFT: ',np.mean((numpy.sin(10 * 2 * numpy.pi * x_recon) - numpy.real(f_recon))**2))
结果：
nfft_inverse
10 iterations: 0.0798988590351581
100 iterations: 0.05136853850272318
1000 iterations: 0.037316315280700896
nfft_gradient_descent
10 iterations: 0.08832834348902704
100 iterations: 0.05901599049633016
1000 iterations: 0.043921864589484
using IFFT: 0.49044932854606377
另一种查看方式是
plt.plot(numpy.sin(10 * 2 * numpy.pi * x_recon),numpy.real(f_recon),'.',label='ifft')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x),numpy.real(nfft_gradient_descent(x,L = 5)),label='gradient descent L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x),numpy.real(nfft_inverse(x,label='nfft_inverse L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x),np.real(f_hat),label='original')
plt.legend()
即使 IFFT 矩阵的条件更好，它也会导致信号重建更差。同样从最后一个图可以看出，有轻微的衰减。可能是由于虚部的系统能量泄漏（我的代码中有错误？？）。只是一个快速测试，乘以 1.3 会得到更好的结果
Bob 已经发布了 an excellent answer，这只是补充一些我希望有启发意义的细节。
首先，比较计算出的频率分量的两个图。请注意 NFFT 的噪声比常规 FFT 的噪声大得多。您正在从噪声样本中估计采样信号的这些频率分量，在一种情况下，样本是规则间隔的，在另一种情况下，它们是随机间隔的。众所周知，常规抽样比随机抽样更有效（有效意味着您需要更少的样本来获得相同数量的信息）。因此，预计随机抽样会产生更多噪声的结果。
我们可以根据 NFFT 估计的频率分量计算“正常”逆 FFT：
f_recon = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k)))
x_recon = numpy.linspace(-.5,.5,N)
我使用ifftshift是因为NFFT定义了k从-N/2到N/2-1，而FFT定义它从0到{{ 1}}。 N-1 交换信号的两半，将第一半变成第二半（ifftshift 从 k 到 N/2 等于 N-1 到 {{1} })。我还在 IFFT 的结果上使用了 -N/2，因为同样的事情适用于时间轴，它将原点从第一个样本移到序列的中间。
请注意 -1 的嘈杂程度。这是由于我们可以对非均匀采样信号做出的 fftshift 的糟糕估计来解释的。还有一个符号错误，当我们比较 f_recon 的两个估计值时，我们已经可以观察到这个错误。这来自于指数中与逆 DFT 具有相同符号的伴随 NFFT，这实际上意味着 f_k 被翻转 w.r.t. x.
如果我们增加随机样本的数量，我们可以获得更好的估计：
f_k
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。