实现了维纳滤波，消除了模糊和噪声图像

Question

我尝试使用以下公式实现维纳滤波器以消除图像模糊：

这里H（k，l）是我的频域模糊核，sigma_v是附加噪声标准差，p_s（k，l）是原始图像的估计功率谱。 我使用以下操作使用15 x 15内核模糊图像：

# set blur size
size_Blur2D = (15,15)
# calculate blur filter coefficients
filter_Blur2D = np.ones(size_Blur2D)
filter_Blur2D /= filter_Blur2D.sum()

将噪声添加为：

# add Gaussian noise to blurred image
noise_Mean        = 0.0
noise_Std         = 2.55
imgBlurred2DNoisy = imgBlurred2D + norm.rvs(noise_Mean, noise_Std, size_img)

然后，我将噪声/模糊（观察到的）图像划分为64个块，每个块的大小为64 x 64，并使用汉宁窗口对其进行窗口处理，以避免DFT瑕疵：

# nSize: window size
# nTaper: taper width
def my_winHanning(nSize, nTaper):

    assert nTaper > 0
    assert nSize > 2 * nTaper

    # create an array to hold the window
    winHann = np.ones(nSize)

    # taper left and right sides
    val = math.pi / nTaper
    for n in range(nTaper):

        # calculate taper
        w = 0.5 + 0.5 * math.cos(n * val)
        # assign right taper
        winHann[nSize - nTaper + n] = w
        # assign left taper
        winHann[nTaper - 1 - n] = w

    return winHann


# expand the blur kernel to power spectrum size
hB, wB             = size_Blur2D
h_expand           = np.zeros((64,64),dtype = 'float64')
h_expand[:hB, :wB] = filter_Blur2D

# cyclicly shift blur kernel to bring its center to (0,0)
hB_half        = hB//2
wB_half        = wB//2
h_expand_shift = np.roll(h_expand, (hB_half, wB_half), (0,0))

# for every 64x64 block of the observe image:
#   window the block with a Hanning window (use a taper length of 6)
#   take the DFT of each block
#   calculate the magnitude-squared of DFT's
#   obtain P_g(k,l) as the average over the blocks (there should be 64 blocks)


def calculate_Pg(img):

    hanning_window = my_winHanning(64,6)
    hanning_window = np.outer(hanning_window, hanning_window)

    # Define each window size
    bl_h            = 64
    bl_w            = 64
    DFT_mag_squares = []
    blocks_DFTs     = []

    for row in np.arange(img.shape[0]  - bl_h + 1, step=bl_h):

        for col in np.arange(img.shape[1]  - bl_w + 1, step=bl_w):

                block = img[row:row+bl_h, col:col+bl_w]
                smoothed_block = block * hanning_window
                block_DFT = ft.fft2(smoothed_block)
                blocks_DFTs.append(block_DFT)
                block_DFT_mag_square = np.abs(block_DFT.real) ** 2 
                DFT_mag_squares.append(block_DFT_mag_square)

    blocks_DFTs             = np.array(blocks_DFTs)            
    DFT_mag_squares         = np.array(DFT_mag_squares)
    print("Number of blocks = {}".format(len(DFT_mag_squares)))
    DFT_mag_squares_sum = 0

    for i in range(0,63):
        DFT_mag_squares_sum += DFT_mag_squares[i]

    P_g = DFT_mag_squares_sum/64 

    return P_g

在找到p_g（块DFT幅值平方的平均值）后，我使用以下约束条件估计p_s：

def calculate_P_s(P_g,blurr_kernel ,noise_std):

    epsilon     = 1e-06
    upper_bound = max(P_g[0,0] - noise_std ** 2 ,epsilon)
    P_s         = []
    H           = ft.fft2(blurr_kernel)

    for i in range(P_g.shape[0]):

        for j in range(P_g.shape[1]):

            first_term = max(P_g[i,j] - noise_std ** 2,epsilon) / max(np.abs(H[i,j]) ** 2 , epsilon)
            P_s.append(min(first_term , upper_bound))
    P_s = np.array(P_s).reshape(64,64)

    return P_s

最后一步是：

# Calculate Phi and its impulse response
factor               = 1
wiener_coeff         = (noise_Std ** 2/P_s_chart)
H                    = ft.fft2(h_expand_shift)
Phi                  = H.conjugate() / ((np.abs(H) ** 2 + factor * wiener_coeff))
Phi_impulseResponse  = (ft.ifft2(Phi)).real
filtered_img         = signal.convolve2d(imgBlurred2DNoisy, Phi_impulseResponse , mode='same', boundary='symm')

但是，获得的p_s和p_g值非常大（例如，最大值可能大于5*10^11），这会导致以下结果：

通过反复尝试，似乎系数*维纳系数应为>；0.001以获得此类结果：

当然，这也不是完美的结果（输出图像由于模糊内核扩展操作而移位）。但至少它表明获得的图像功率谱太大，无法用该滤波器实现。我已经在互联网上查看了维纳滤波器的实现，但是没有找到一个像我正在尝试的那样详细的实现。这些公式或我的代码有什么问题吗