加性高斯白噪声及维纳滤波的基本原理与Python实现

2019-05-25 本文已影响0人 iwuqing

1. 加性高斯白噪声

1.1 特点

加性高斯白噪声属于白噪声的一种，有如下两个特点：

功率谱密度为常数，即白
高斯则为其幅度的概率密度函数服从高斯分布

1.2 高斯分布的一维概率密度函数

$P(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$

$\sigma$ 为标准差
$\mu$ 为均值

1.3 Python内置库说明

random.gauss(mu, sigma)其值即服从高斯分布，若想要是实现加性高斯白噪声，循环作加即可

1.4 Code

import numpy as np
import random

def gaussian_white_noise(intput_signal, mu, sigma):
    '''
    加性高斯白噪声（适用于灰度图）
    :param intput_signal: 输入图像
    :param mu: 均值
    :param sigma: 标准差
    :return:
    '''

    intput_signal_cp = np.copy(intput_signal)  # 输入图像的副本

    m, n = intput_signal_cp.shape   # 输入图像尺寸（行、列）

    # 添加高斯白噪声
    for i in range(m):
        for j in range(n):
            intput_signal_cp[i, j] = intput_signal_cp[i, j] + random.gauss(mu, sigma)


    return intput_signal_cp

1.5 测试对比

均值为0，标准差为10时的加行高斯白噪声前后对比图

2. 逆滤波

实际上逆滤波是维纳滤波的一种理想情况，当不存在加性噪声时，维纳滤波与逆滤波等同。

2.1 基本原理

在时域内有
$y(n) = x(n)\bigotimes h(n)$

$x(n)$ 为输入信号
$h(n)$ 为系统函数
$y(n)$ 为输出信号

根据时域卷积定理，我们知道时域卷积等于频域乘积

$Y(\omega) = X(\omega) \times H(\omega)$

$Y, X, H$ 分别为输出信号，输入信号，系统函数的傅里叶变换

则有
$X(\omega) = \frac{Y(\omega)}{H(\omega)}$
这意味着，当我们已知系统函数时，我们可以很简单的完成滤波。

2.2 Code

import numpy as np


def inverse_filtering(input_singal, h):
    '''
    逆滤波
    :param input_singal: 输入信号
    :param h: 退化函数（时域）
    :return: 滤波后的信号（幅值）
    '''
    output_signal = [] #输出信号
    output_signal_fft = [] #输出信号的傅里叶变换

    h_fft = np.fft.fft2(h) # 退化函数的傅里叶变换

    input_singal_fft = np.fft.fft2(input_singal)  # 输入信号的傅里叶变换

    output_signal_fft = input_singal_fft / h_fft

    output_signal = np.abs(np.fft.ifft2(output_signal_fft))

    return output_signal

3 维纳滤波

3.1 直观认识

理解了逆滤波的基本过程之后，实际上维纳滤波就不是太大问题了。实际上，逆滤波对于绝大多数情况滤波效果都不好，因为逆滤波是通过傅里叶变换将信号由时域转换到频域，再根据时域卷积定理，在频域作除法。对于乘性干扰这当然是没问题的，甚至是完美的。而如果存在加性噪声，例如：加性高斯白噪声。逆滤波效果就不好了，某些情况下几乎无法完成滤波情况。

3.2 数学解释

输入信号经过系统函数后
时域上
$y(n) = x(n)\bigotimes h(n)$
频域上
$Y(\omega) = X(\omega) \times H(\omega)$
若存在加性噪声则为
时域上
$y(n) = g(n) + d(n)$

其中 $g(n) = x(n)\bigotimes h(n)$

频域上
$Y(\omega) = G(\omega) + C(\omega)$

其中 $G(\omega) = X(\omega) \times H(\omega)$

则
$\frac {Y(\omega)}{H(\omega)} = \hat{X}(\omega)=X(\omega)+\frac{C(\omega)}{H(\omega)}$

于是，从上面对输入信号的估计表达式可以看出，多出了一项加性噪声的傅里叶变换与系统函数的比值。尤其当 $H(\omega)$ 相对于 $C(\omega)$ 很小时，滤波后的信号差距十分严重。

3.3 维纳滤波表达式

$\hat{X}(\omega)=\left[\frac{1}{H(\omega)} \frac{|H(\omega)|^{2}}{|H(\omega)|^{2}+K}\right] G(\omega)$

$K = \frac{|C(\omega)|^2}{|H(\omega)|^2}$

3.4 Code

import numpy as np

def wiener_filtering(input_signal, h, K):
    '''
    维纳滤波
    :param input_signal: 输入信号
    :param h: 退化函数（时域）
    :param K: 参数K
    :return: 维纳滤波后的信号（幅值）
    '''
    output_signal = []  # 输出信号

    output_signal_fft = [] # 输出信号的傅里叶变换

    input_signal_cp = np.copy(input_signal) # 输入信号的副本

    input_signal_cp_fft = np.fft.fft2(input_signal_cp)  # 输入信号的傅里叶变换

    h_fft = np.fft.fft2(h) # 退化函数的傅里叶变换

    h_abs_square = np.abs(h_fft)**2 # 退化函数模值的平方

    # 维纳滤波
    output_signal_fft = np.conj(h_fft) / (h_abs_square + K)

    output_signal = np.abs(np.fft.ifft2(output_signal_fft * input_signal_cp_fft)) # 输出信号傅里叶反变换

    return output_signal

4 维纳滤波与逆滤波对比

camera

coins

5 GitHub

click me!