示例 python nfft 傅里叶变换 - 信号重建归一化问题
Example python nfft fourier transform - Issues with signal reconstruction normalization
我为两者写了一个完整的工作示例 nfft, and 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, f, 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, .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,-499,...-1]
f_k = scipy.fft.fft(f)
#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 scipy.fft
f_recon = scipy.fft.ifft(f_k)
#Plot original vs reconstructed
plt.figure()
plt.title('scipy.fft')
plt.scatter(x, f, label='f(x)')
plt.scatter(x, f_recon, label='f_recon(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, N, L = 1)
plt.hist(f_hat - numpy.real(h_hat), bins=30, label='1 iteration')
h_hat = nfft_inverse(x, f, N, L = 10)
plt.hist(f_hat - numpy.real(h_hat), bins=30, label='10 iterations')
h_hat = nfft_inverse(x, f, N, L = 1000)
plt.hist(f_hat - numpy.real(h_hat), bins=30, 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, y, N, 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, r, 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
Is it basically true that you found that the there is not a true inverse to the algorithm then? I cannot obtain my original points? x != nfft(nfft_adjoint(x))
请检查参考文献的第 2.3 节 paper
数值比较
提到了另一种可能性，即不是在x点重建f，你可以使用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, f, 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, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_inverse(x, f, len(x), 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, f, len(x), L = 10)
print('10 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 1000)
print('1000 iterations: ', 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, f, len(x), L = 5)), '.', label='gradient descent L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x), numpy.real(nfft_inverse(x, f, len(x), L = 5)), '.', 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 已经发布 ，这只是为了补充一些细节，希望对您有所启发。
首先，比较计算出的频率分量的两个图。请注意 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 到N-1。 ifftshift 交换信号的两半，将第一部分变成第二部分（k 从 N/2 到 N-1 等于 -N/2 到 -1).我还在 IFFT 的结果上使用了 fftshift，因为同样的事情适用于时间轴，它将原点从第一个样本移动到序列的中间。
请注意 f_recon 的噪音有多大。这是由于我们可以用非均匀采样信号做出的 f_k 的不良估计来解释的。还有一个符号错误，当我们比较 f_k 的两个估计值时，我们已经可以观察到这个错误。这来自与逆 DFT 在指数中具有相同符号的伴随 NFFT，这实际上意味着 f_recon 被翻转 w.r.t。 x.
如果我们增加随机样本的数量，我们可以获得更好的估计：
import numpy
import nfft
import matplotlib.pyplot as plt
#Define signal:
N = 1024 * 16 # power of two for speed
x = -0.5 + numpy.random.rand(N)
f = numpy.sin(10 * 2 * numpy.pi * x) + .1 * numpy.random.randn(N) # Add some 'y' randomness to the sample
#prepare wavenumbers for transform:
k = - N // 2 + numpy.arange(N)
N2 = 1024
f_k = nfft.nfft_adjoint(x, f, N2, truncated=False)
#Reconstruct the original signal with nfft
# (note the minus sign to flip the signal, in reality we should flip x)
f_recon = - numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k))) / (N / N2)
x_recon = numpy.linspace(-.5, .5, N2, endpoint=False)
#Plot original vs reconstructed
plt.figure()
plt.title('nfft')
plt.scatter(x[:N2], f[:N2], label='f(x)') # don't plot all samples, there's too many
plt.scatter(x_recon, f_recon, label='f_recon(x)')
plt.legend()
plt.show()
我为两者写了一个完整的工作示例 nfft, and 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, f, 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, .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,-499,...-1]
f_k = scipy.fft.fft(f)
#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 scipy.fft
f_recon = scipy.fft.ifft(f_k)
#Plot original vs reconstructed
plt.figure()
plt.title('scipy.fft')
plt.scatter(x, f, label='f(x)')
plt.scatter(x, f_recon, label='f_recon(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, N, L = 1)
plt.hist(f_hat - numpy.real(h_hat), bins=30, label='1 iteration')
h_hat = nfft_inverse(x, f, N, L = 10)
plt.hist(f_hat - numpy.real(h_hat), bins=30, label='10 iterations')
h_hat = nfft_inverse(x, f, N, L = 1000)
plt.hist(f_hat - numpy.real(h_hat), bins=30, 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, y, N, 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, r, 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
Is it basically true that you found that the there is not a true inverse to the algorithm then? I cannot obtain my original points? x != nfft(nfft_adjoint(x))
请检查参考文献的第 2.3 节 paper
数值比较
x点重建f，你可以使用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, f, 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, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_inverse(x, f, len(x), 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, f, len(x), L = 10)
print('10 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 100)
print('100 iterations: ', np.mean((numpy.sin(10 * 2 * numpy.pi * x) - numpy.real(h_hat))**2))
h_hat = nfft_gradient_descent(x, f, len(x), L = 1000)
print('1000 iterations: ', 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, f, len(x), L = 5)), '.', label='gradient descent L=5')
plt.plot(numpy.sin(10 * 2 * numpy.pi * x), numpy.real(nfft_inverse(x, f, len(x), L = 5)), '.', 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 已经发布
首先，比较计算出的频率分量的两个图。请注意 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 到N-1。 ifftshift 交换信号的两半，将第一部分变成第二部分（k 从 N/2 到 N-1 等于 -N/2 到 -1).我还在 IFFT 的结果上使用了 fftshift，因为同样的事情适用于时间轴，它将原点从第一个样本移动到序列的中间。
请注意 f_recon 的噪音有多大。这是由于我们可以用非均匀采样信号做出的 f_k 的不良估计来解释的。还有一个符号错误，当我们比较 f_k 的两个估计值时，我们已经可以观察到这个错误。这来自与逆 DFT 在指数中具有相同符号的伴随 NFFT，这实际上意味着 f_recon 被翻转 w.r.t。 x.
如果我们增加随机样本的数量，我们可以获得更好的估计：
import numpy
import nfft
import matplotlib.pyplot as plt
#Define signal:
N = 1024 * 16 # power of two for speed
x = -0.5 + numpy.random.rand(N)
f = numpy.sin(10 * 2 * numpy.pi * x) + .1 * numpy.random.randn(N) # Add some 'y' randomness to the sample
#prepare wavenumbers for transform:
k = - N // 2 + numpy.arange(N)
N2 = 1024
f_k = nfft.nfft_adjoint(x, f, N2, truncated=False)
#Reconstruct the original signal with nfft
# (note the minus sign to flip the signal, in reality we should flip x)
f_recon = - numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(f_k))) / (N / N2)
x_recon = numpy.linspace(-.5, .5, N2, endpoint=False)
#Plot original vs reconstructed
plt.figure()
plt.title('nfft')
plt.scatter(x[:N2], f[:N2], label='f(x)') # don't plot all samples, there's too many
plt.scatter(x_recon, f_recon, label='f_recon(x)')
plt.legend()
plt.show()