小波变换的前因后果（三）

小波去噪是建立在DWT的基础上的，需要进行小波分解、再重构。接上一篇。

小波分析即用Mallat塔式算法对信号进行降阶分解。该算法在每尺度下将信号分解成近似分量与细节分量。近似分量表示信号的高尺度，即低频信息；细节分量表示信号的低尺度，即高频信息。对含有噪声的信号，噪声分量的主要能量集中在小波分解的细节分量中。

 二、小波去噪

 1、概念

通常情况下， 我们在从设备上采集到的信号都是具有一定的噪声的，大多数情况下，可认为这种噪声为高斯白噪声。被噪声污染的信号=干净的信号+噪声。

为什么要使用阈值：由于信号在空间上(或者时间域)是有一定连续性的，因此在小波域，有效信号所产生的小波系数其模值往往较大；而高斯白噪声在空间上(或者时间域)是没有连续性的，因此噪声经过小波变换，在小波阈仍然表现为很强的随机性，通常仍认为是高斯白噪的。 那么就得到这样一个结论：在小波域，有效信号对应的系数很大，而噪声对应的系数很小。 刚刚已经说了，噪声在小波域对应的系数仍满足高斯白噪分布。如果在小波域，噪声的小波系数对应的方差为sigma，那么根据高斯分布的特性，绝大部分(99.99%)噪声系数都位于[-3*sigma，3*sigma]区间内(切比雪夫不等式， 3sigma准则)。因此，只要将区间[-3*sigma，3*sigma]内的系数置零(这就是常用的硬阈值函数的作用)，就能最大程度抑制噪声的，同时只是稍微损伤有效信号。将经过阈值处理后的小波系数重构，就可以得到去噪后的信号。 常用的软阈值函数，是为了解决硬阈值函数“一刀切”导致的影响(模小于3*sigma的小波系数全部切除，大于3*sigma全部保留，势必会在小波域产生突变，导致去噪后结果产生局部的抖动，类似于傅立叶变换中频域的阶跃会在时域产生拖尾)。软阈值函数将模小于3*sigma的小波系数全部置零，而将模大于3*sigma的做一个比较特殊的处理，大于3*sigma的小波系数统一减去3*sigma，小于-3*sigma的小波系数统一加3*sigma。经过软阈值函数的作用，小波系数在小波域就比较光滑了，因此用软阈值去噪得到的图象看起来很平滑，类似于冬天通过窗户看外面一样，像有层雾罩在图像上似的。

比较硬阈值函数去噪和软阈值函数去噪：硬阈值函数去噪所得到的峰值信噪比(PSNR)较高，但是有局部抖动的现象；软阈值函数去噪所得到的PSNR不如硬阈值函数去噪，但是结果看起来很平滑，原因就是软阈值函数对小波系数进行了较大的 “社会主义改造”，小波系数改变很大。因此各种各样的阈值函数就出现了，其目的我认为就是要使大的系数保留，小的系数被剔出，而且在小波域系数过渡要平滑。

如何估计小波域噪声方差sigma的估计，这个很简单：把信号做小波变换，在每一个子带利用robust estimator估计就可以(可能高频带和低频带的方差不同)。 robust estimator就是将子带内的小波系数模按大小排列，然后取最中间那个，然后把最中间这个除以0.6745就得到噪声在某个子带内的方差sigma。利用这个sigma，然后选种阈值函数，就可以去去噪了，在matlab有实现api可使用。

小波阈值去噪过程如下图

在小波分析中经常用到近似和细节，近似表示信号的高尺度，即低频信息；细节表示信号的低尺度，即高频信息。对含有噪声的信号，噪声分量的主要能量集中在小波解的细节分量中。 

2、原理

小波阈值去噪的实质为抑制信号中无用部分、增强有用部分的过程。小波阈值去噪过程为：(1)分解过程，即选定一种小波对信号进行n层小波分解；(2)阈值处理过程，即对分解的各层系数进行阈值处理，获得估计小波系数；(3)重构过程，据去噪后的小波系数进行小波重构，获得去噪后的信号。

小波阈值去噪过程

小波分解重构过程

小波分解：X->ca3,cd3,cd2,cd1；小波重构：ca3,cd3,cd2,cd1->X。其中ca为低频信息、近似分量，cd为高频、细节分量。

3、影响降噪效果的因素

（1）小波基的选择

在对信号进行小波分解时需要选择合适的小波基，由于没有任何一种小波基可以对不同类型的信号达到最优的分解效果，因此，如何选择小波基成为小波分解的一个重点。针对现实中的信号，小波基的选择一般要考虑以下几个因素：支撑长度、对称性、消失矩、正则性、相似性。针对一维信号，例如语音信号，通常选择dB小波和sym小波。

（2）分解层数的选择

在对信号进行小波分解时，分解的层数取得越大，则噪声和信号表现的不同特性越明显，越有利于二者的分离，但是分解的层数越大，经过重构的信号失真也会越大，在一定程度上会对信号去噪的效果产生较差的影响。因此，如何选择分解层数以解决信噪分离效果和重构信号失真之间的矛盾呢？

小波分解的频段范围与采样频率有关。若进行N层分解，则各个频段范围为：

假设原始信号X的采样频率为1000Hz，则信号的最大频率为500，对该信号做3层小波分解，则各个频段范围如下图所示。

通常小波分解的频段范围与采样频率有关。若N层分解，则各个频段大小为Fs/2/2^N 。例如：一个原始信号，经历的时间长度为2秒，采样了2000个点，那么做除法，可得出采样频率为1000hz，由采样定理（做除法）得该信号的最大频率为500hz，那么对该信号做3层的DWT，一阶细节的频段为250－500hz，一阶逼近的频段为小于250hz，二阶细节的频段为125－250hz，逼近的频段为小于125hz，三阶细节的频段约为62.5－125hz，逼近的频段为小于62.5hz。对于更多阶的分解也是以此类推的。

 （3）具体的分解算法（Mallat塔式算法）

多分辨率分析为正交小波基的构造提供了一种简单方法，同时它还是正交小波变换的快速算法(即 Mallat 算法)基本理论基础。Mallat 算法是由 S.Mallat于1989年提出的，该算法在小波分析中的作用相当于快速傅立叶变换在傅里叶分析中的作用。Mallat 算法由小波滤波器 H、G、和 h、g 对测量的信号进行分解和重构。它的分解算法可以表述为：

分解算法可以用图解的形式表示为：

式中，h，g为时域中的小波重构滤波器，j为分解的层数，若分解的深度(分解的最高层)为J，则j=J−1,J−2,…,1,0，Aj为信号f(t)在第j层的低频部分(近似部分)的小波系数；Dj为信号f(t)在第 j 层的高频部分(细节部分)的小波系数。 

 这里提供了一维信号下的小波分解代码

"""
将数据序列进行小波分解，每一层分解的结果是上次分解得到的低频信号再分解成低频和高频两个部分。
如此进过N层分解后源信号X被分解为：X = D1 + D2 + … + DN + AN 
其中D1,D2,…,DN分别为第一层、第二层到等N层分解得到的高频信号，AN为第N层分解得到的低频信号。
"""
import numpy as np
import matplotlib.pyplot as plt

import pywt
import pywt.data


ecg = pywt.data.ecg()
print(len(ecg))
data1 = np.concatenate((np.arange(1, 400),
                        np.arange(398, 600),
                        np.arange(601, 1024)))
x = np.linspace(0.082, 2.128, num=1024)[::-1]
data2 = np.sin(40 * np.log(x)) * np.sign((np.log(x)))

mode = pywt.Modes.smooth


def plot_signal_decomp(data, w, title):
    """Decompose and plot a signal S.
    S = An + Dn + Dn-1 + ... + D1
    """
    w = pywt.Wavelet(w)#选取小波函数
    a = data
    ca = []#近似分量
    cd = []#细节分量
    for i in range(5):
        (a, d) = pywt.dwt(a, w, mode) #进行5阶离散小波变换
        ca.append(a)
        cd.append(d)

    rec_a = []
    rec_d = []

    for i, coeff in enumerate(ca):
        coeff_list = [coeff, None] + [None] * i
        rec_a.append(pywt.waverec(coeff_list, w)) #重构

    for i, coeff in enumerate(cd):
        coeff_list = [None, coeff] + [None] * i
        if i == 3:
            print(len(coeff))
            print(len(coeff_list))
        rec_d.append(pywt.waverec(coeff_list, w))

    fig = plt.figure()
    ax_main = fig.add_subplot(len(rec_a) + 1, 1, 1)
    ax_main.set_title(title)
    ax_main.plot(data)
    ax_main.set_xlim(0, len(data) - 1)

    for i, y in enumerate(rec_a):
        ax = fig.add_subplot(len(rec_a) + 1, 2, 3 + i * 2)
        ax.plot(y, 'r')
        ax.set_xlim(0, len(y) - 1)
        ax.set_ylabel("A%d" % (i + 1))

    for i, y in enumerate(rec_d):
        ax = fig.add_subplot(len(rec_d) + 1, 2, 4 + i * 2)
        ax.plot(y, 'g')
        ax.set_xlim(0, len(y) - 1)
        ax.set_ylabel("D%d" % (i + 1))


# plot_signal_decomp(data1, 'coif5', "DWT: Signal irregularity")
# plot_signal_decomp(data2, 'sym5',
#                   "DWT: Frequency and phase change - Symmlets5")
plot_signal_decomp(ecg, 'sym5', "DWT: Ecg sample - Symmlets5")

plt.show()
"""
将数据序列进行小波分解，每一层分解的结果是上次分解得到的低频信号再分解成低频和高频两个部分。
如此进过N层分解后源信号X被分解为：X = D1 + D2 + … + DN + AN 
其中D1,D2,…,DN分别为第一层、第二层到等N层分解得到的高频信号，AN为第N层分解得到的低频信号。
"""

二维信号的小波分解代码如下

# 采用小波阈值方法对预先加了标准差为15的高斯噪声图像(256*256)进行小波阈值去噪，
# 其中各层的滤波阈值均为同一值，采用的小波基为sym4小波，分解层数为3，阈值函数为软阈值函数

"""
原理：
小波阈值去噪的是指为抑制信号中无用部分、增强有用部分的过程，小波阈值去噪过程为：
1、分解过程，即选定一种小波对信号进行n层小波分解
2、阈值处理过程，即对分解的隔层系数进行阈值处理，获得估计小波系数
3、重构过程，据去噪后的小波系数进行小波重构，获得去噪后的信号

影响效果的因素：分解层数、阈值、小波基的选择、阈值函数的选择
"""

import numpy as np
import pywt
from cv2 import COLOR_BGR2GRAY, COLOR_RGB2GRAY, cv2, denoise_TVL1
from PIL import Image

# img = cv2.imread("lenage15.bmp",COLOR_BGR2GRAY)
img = cv2.imread("lenage15.bmp",0)
print(img.shape)
w = 'sym4' #定义小波基的类型
l = 3 # 变换层次为3
coeffs = pywt.wavedec2(data=img, wavelet=w, level=l) # 对图像进行小波分解
"""
coeffs shape：
4*38*38
"""
threshold = 0.04

list_coeffs = []
for i in range(1, len(coeffs)):
    list_coeffs_ = list(coeffs[i]) # 38*38
    list_coeffs.append(list_coeffs_)
# print(list_coeffs) # list_coeffs.shape:4*38*38

for r1 in range(len(list_coeffs)): # 4
    for r2 in range(len(list_coeffs[r1])): # 38
        # 对噪声滤波（软阈值）
        list_coeffs[r1][r2] = pywt.threshold(list_coeffs[r1][r2], threshold*np.max(list_coeffs[r1][r2]))

rec_coeffs = [] # 重构系数
rec_coeffs.append(coeffs[0]) #将原图像的低尺度系数保留进来

for j in range(len(list_coeffs)):
    rec_coeffs_ = tuple(list_coeffs[j])
    rec_coeffs.append(rec_coeffs_)


denoised_img = pywt.waverec2(rec_coeffs, 'sym4')
denoised_img = Image.fromarray(np.uint8(denoised_img))
denoised_img.save("result.bmp")
# 总结：这里能够看到小波去噪对去除高斯噪声的效果还是可以的，但是效果不是很明显。
# 大家可以从阈值的选择以及调整滤波的阈值，甚至将各层的阈值设置为不同的值再去试试。
# 再看看有关的论文关于提升小波阈值效果的实现方法。# 本文只是举个简单的例子来说明小波阈值的处理过程