逆向过滤可消除卷积爆炸

Question

本质上，我正在写一篇论文，其中我想要一个图来显示任意曲线与高斯卷积的效果。然后，我想展示您可以通过将卷积形式的 FFT 除以高斯的 FFT（逆滤波）来对其进行反卷积。然而这样做却会引起爆炸，我也不知道为什么。

在说之前，我知道逆滤波是一种糟糕的方法，而且对噪声非常敏感。然而，这就是这个数字的要点。基本上表明它有效，但一旦你离开非理想条件它就会爆炸

首先，我设置了基本导入并定义了它们所在的函数/域。

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from scipy.fft import *


xf = 5
n = 100
xrange = np.linspace(0,xf,num=n)

a = 10
xc = 2
w = 0.1

def func(x):
    bpoint = 2.25
    if x <= bpoint:
        return a/(np.pi*w*(1+((x-xc)**2/w)))
    else:
        return a/(np.pi*w*(1+((x-xc)**2/w))) + 9*np.sin(1.5*(x-bpoint)) + 5*(x - bpoint)
    
def gauss(x,mx,o,I):
    return I*np.exp(-1/2*(((x-mx)/(o/6))**2))

def conved(x,mx,o,I):
    return gauss(x,mx,o,I)*func(x)

然后，我将函数 f(x) 与要与之卷积的高斯 g(x) 一起绘制。

f = np.zeros((len(xrange),))
for i in range(len(xrange)):
    f[i] = func(xrange[i])

o = 2.5

g = gauss(xrange,2.5,o,5)

normal = np.sqrt(2*np.pi)*o/6*5

plt.figure(0)
plt.plot(xrange,f)
plt.plot(xrange,g)
plt.legend(['f(x)','g(x)'])

从那里，我使用积分公式（经过归一化）计算卷积，并将它们并排绘制。

fconvg = np.zeros((len(xrange),))
for i in range(len(xrange)):
    mx = xrange[i]
    avg = sp.integrate.quad(conved,mx-o/2,mx+o/2,args=(mx,o,5))
    fconvg[i] = avg[0]/normal

plt.figure(1)
plt.plot(xrange,f)
plt.plot(xrange,fconvg)
plt.legend(['f(x)','f(x)*g(x)'])

现在我知道将每个函数的 FFT 相乘可以对它们进行卷积。然而，FFT 处理边界条件的方式并不准确，如下图所示。

F = fft(f)
G = fft(g)

fconvg_fft = fftshift(ifft(F*G))*xf/(n-1)/normal

plt.figure(2)
plt.plot(xrange,f)
plt.plot(xrange,fconvg)
plt.plot(xrange,fconvg_fft)
plt.legend(['f(x)','Integral','FFT'])

由于它应用了周期性 BC，因此卷积的边界与使用积分计算的真实卷积不同

这部分是我遇到问题的地方。基本上，我采用两种不同卷积形式的 FFT，然后将它们除以高斯的 FFT。然而，在这两种情况下，当我进行 IFFT 时，结果都会膨胀到无穷大。我不知道发生了什么事。在我的原始代码中，这对于使用 FFT 计算卷积的情况非常有效。但在这个例子中，这两种情况都不起作用

H_integral = fft(fconvg)
F_integral = H_integral/G
f_integral = ifft(F_integral)

H_fft = fft(fconvg_fft)
F_fft = H_fft/G
f_fft = ifft(F_fft)

plt.figure(3)
plt.plot(xrange,f)
plt.plot(xrange,f_integral)
plt.legend(['f(x)','Integral'])

plt.figure(4)
plt.plot(xrange,f)
plt.plot(xrange,f_fft)
plt.legend(['f(x)','FFT'])
plt.show()

积分计算反卷积

通过 FFT 计算进行反卷积

我尝试过使用步长进行缩放，并使用 fftshift/ifftshift 进行操作，但没有任何效果。我假设我只是不明白 FFT 是如何工作的，而且它开始变得非常令人沮丧。如果有人能找出我做错了什么，我真的很感激

编辑：由于有人抱怨代码，所以这就是一个块

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from scipy.fft import *

xf = 5
n = 100
xrange = np.linspace(0,xf,num=n)

a = 10
xc = 2
w = 0.1

def func(x):
    bpoint = 2.25
    if x <= bpoint:
        return a/(np.pi*w*(1+((x-xc)**2/w)))
    else:
        return a/(np.pi*w*(1+((x-xc)**2/w))) + 9*np.sin(1.5*(x-bpoint)) + 5*(x - bpoint)
    
def gauss(x,mx,o,I):
    return I*np.exp(-1/2*(((x-mx)/(o/6))**2))

def conved(x,mx,o,I):
    return gauss(x,mx,o,I)*func(x)

f = np.zeros((len(xrange),))
for i in range(len(xrange)):
    f[i] = func(xrange[i])

o = 2.5

g = gauss(xrange,2.5,o,5)

normal = np.sqrt(2*np.pi)*o/6*5

fconvg = np.zeros((len(xrange),))
for i in range(len(xrange)):
    mx = xrange[i]
    avg = sp.integrate.quad(conved,mx-o/2,mx+o/2,args=(mx,o,5))
    fconvg[i] = avg[0]/normal

F = fft(f)
G = fft(g)

fconvg_fft = fftshift(ifft(F*G))*xf/(n-1)/normal

H_integral = fft(fconvg)
F_integral = H_integral/G
f_integral = ifft(F_integral)

H_fft = fft(fconvg_fft)
F_fft = H_fft/G
f_fft = ifft(F_fft)

>>> F_integral = H_integral/G
<ipython-input-39-34c973d32c2b>:2: RuntimeWarning:

divide by zero encountered in divide

<ipython-input-39-34c973d32c2b>:2: RuntimeWarning:

invalid value encountered in divide
>>> G[50]
-0j
>>> G[50] == 0
True

这就是为什么它“炸到无穷大”。如果您要除以这些值，则需要处理 FFT 结果为零或非常接近零的情况。

Answer 1