确定性傅里叶反褶积

Question

我用numpy做傅里叶反褶积有点麻烦。我目前正在尝试用一个3高斯函数的测试用例，这样我就知道在每一端都会发生什么。在

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

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

#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. 这显然是一个规模问题。我唯一能想到这一点的地方是我应用天真约束的时候。我应该在这一步中重新规范化吗，知道fourier空间上的1/sqrt（N）*积分等于1？在
2. 看起来恢复的高斯曲线的位置和曲线两边的一半都弄乱了。这是因为Fourier空间的除法吗？我如何恢复原来的位置（或者我一开始就完全错了）

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

干杯， 基文

编辑：我应该补充我没有问题恢复图像使用杂乱的+一些小的编辑。这一定意味着我的方法有点不对劲吧？在

Answer 1

1）正如您正确理解的，缩放要求与离散傅里叶变换有关。最好的方法是计算两个均匀单位信号的反褶积。它们的DFT是n0 0 0….，其中n是DFT的点数。因此比值r碜f为10 0 0 0，由np.fft.ifft()计算的后向fft为1/n1/n1/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()可用于此目的：

^{pr2}$

编辑：为了研究转换的原因，让我们关注一下反褶积核是一个非常窄的高斯分布，几乎与狄拉克分布相对应。你的输入信号是一条以0为中心的高斯曲线，被采样的帧在-5和5之间。同样，反褶积核是以零为中心的狄拉克。因此，反褶积信号必须与输入信号相同：以零为中心的高斯曲线。然而，在FFTW中实现的DFT，因此np.fft.fft()is computed as that of a frame starting at 0 and ending at 10，在点10j/n取样，其中j在[0..n-1]中，傅里叶空间中的频率为k/10，其中k在[0..n/2，-n/2+1..-1]。因此，这个DFT将信号视为以5为中心的高斯信号，而反褶积内核则视为以5为中心的狄拉克。函数f（t）与以t0为中心的狄拉克三角洲（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：这种滤波器可以设计成在给定卷积信号的信噪比的情况下使反褶积信号的均方误差最小。然而，它需要一些关于信号和噪声的功率谱密度的假设。在