确定性傅里叶反卷积

Question

我在使用 numpy 进行傅里叶反卷积时遇到了一些麻烦。我目前正在尝试使用 3 个高斯的测试用例，所以我确切地知道每一端会发生什么。

我要恢复的是输入信号，给出了滤波器和输出的确切形式。

在这里，我使用了一个简单的约束来抑制高频端，将其设置为零（因为信号也是傅立叶空间中的高斯信号）。由于这个限制，我希望通过一点点振铃来恢复我的原始输入。

#Dummy Case for Gaussian convolve with Gaussian

N = 128
x = np.arange(-5, 5, 10./(2 * N))
epsilon = 1e-18

def gaus(x,sigma):
    return 1./np.sqrt(2*np.pi)/sigma * np.exp(-(x * x)/(2 * sigma**2))

y_g = gaus(x,0.3)     #output gaussian blurred signal
y_b = gaus(x,0.1)     #gaussian blur filter
y_i = gaus(x,np.sqrt(0.3**2 - 0.1**2))   #og gaussian input

f_yg = np.fft.fft(y_g)    #fft the blur
f_yb = np.fft.fft(y_b)    #fft the filter
f_yi = np.fft.fft(y_i)

r_f = (np.fft.fftshift(f_yg)+epsilon)/(np.fft.fftshift(f_yb)+epsilon)      #deconvolve by division in fourier space
r_f[np.abs(x)>0.5] = 0    #naive constraint to remove the artifacts by knowing final form is gaussian

r_f = np.fft.ifftshift(r_f)

r_if = np.fft.ifft(r_f)
y_gf = np.fft.ifft(f_yg)
y_bf = np.fft.ifft(f_yb)
y_if = np.fft.ifft(f_yi)

plt.plot(x,y_if, label='fft true input')
plt.plot(x,r_if, label='fft recv. input')
plt.legend(framealpha=0.)
plt.show()

这里的橙色是使用输出和模糊的反卷积恢复的输入信号。

我有几个问题：

1. 显然存在缩放问题。我认为这可能会出现的唯一领域是当我应用天真约束时。如果知道傅立叶空间上的 1/sqrt(N)*integral 等于 1，我是否应该在这一步重新归一化？
2. 看起来恢复的高斯的位置与曲线两侧的一半曲线混淆了。这是由于傅立叶空间的划分吗？我如何恢复原来的位置（或者我一开始就做错了）

我附上了用于生成两条曲线的脚本，即物理空间中的原始输入和恢复输入。

干杯， 凯文

编辑：我应该补充一点，使用 scipy.deconvolve + 一些小编辑恢复图像没有问题。这一定意味着我这里的方法有些错误？

Answer 1

1 ) 正如您正确理解的那样，缩放的要求与离散傅里叶变换有关。获得它的最好方法是计算两个统一单位信号的反卷积。它们的 DFT 为 n 0 0 0 ....，其中 n 是 DFT 的点数。因此比率 r_f 是 1 0 0 0 0 并且由np.fft.ifft() 计算的反向 fft 是 1/n 1/n 1/n ... 反卷积产生的正确信号应该是 1/T 1/T 1/T ...，其中 T=10。是帧的长度。

因此，执行反卷积的正确缩放比例是n/T= len(r_f)/10.

r_if=r_if*len(r_if)/10.

2) 去卷积信号被平移半个周期。这是因为高斯核位于帧的中间。只需将内核移动半个周期即可解决问题。函数np.fft.fftshift()可以应用于此：

f_yb = np.fft.fft(np.fft.fftshift(y_b))    #fft the filter

编辑：为了调查翻译的原因，让我们关注反卷积核是非常窄的高斯分布的情况，几乎对应于狄拉克分布。您的输入信号是一条高斯曲线，以零为中心，帧在 -5 和 5 之间进行采样。类似地，反卷积内核是一个以零为中心的狄拉克。因此，去卷积信号必须与输入信号相同：以零为中心的高斯曲线。尽管如此，在 FFTW 中实现的 DFT 以及因此在 np.fft.fft() is computed as that of a frame starting at 0 and ending at 10 中实现的 DFT，在点 10j/n 处采样，其中 j 在 [0..n-1] 中，傅立叶空间中的频率为 k/10，其中 k 在 [ 0..n/2,-n/2+1..-1]。因此，此 DFT 将您的信号视为以 5 为中心的高斯，将反卷积核视为以 5 为中心的 Dirac。函数 f(t) 与以 t0 为中心的 Dirac delta(t-t0) 的卷积就是翻译函数 f(t-t0)。因此，np.fft.fft() 计算的反卷积结果是输入信号平移了半个周期。由于输入信号在 [-5,5] 帧中以 0 为中心，因此np.fft.fft() 计算的输出信号以 -5 为中心（或由于周期性而等效为 5）。 移动内核解决了我们将帧视为 [-5 5] 和 np.fft.ifft() 将其视为 [0 10] 处理之间的不匹配。

滤波器通常旨在减少高频噪声的影响。因此，去卷积会导致高频噪声的潜在放大。像您一样筛选频率是一个潜在的解决方案。请注意，它完全等同于使用特定滤波器对信号进行卷积！

在断层重建的范围内，滤波后的反投影算法需要应用斜坡滤波器，这会极大地放大高频噪声。 Here 提出了一个Wiener filter：在给定卷积信号的 SNR 的情况下，可以设计这种滤波器以最小化去卷积信号的均方误差。然而，它需要一些关于信号和噪声的功率谱密度的假设。