1、前言
在现实生活和工做中，噪声无处不在，在许多领域中，如天文、医学图像和计算机视觉方面收集到的数据经常是含有噪声的。噪声可能来自获取数据的过程，也可能来自环境影响。因为种种缘由,总会存在噪声，噪声的存在每每会掩盖信号自己所要表现的信息，因此在实际的信号处理中，经常须要对信号进行预处理，而预处理最主要的一个步骤就是降噪。html
小波分析是近年来发展起来的一种新的信号处理工具，这种方法源于傅立叶分析，小波（wavelet），即小区域的波，仅仅在很是有限的一段区间有非零值，而不是像正弦波和余弦波那样无始无终。小波能够沿时间轴先后平移，也可按比例伸展和压缩以获取低频和高频小波，构造好的小波函数能够用于滤波或压缩信号，从而能够提取出已含噪声信号中的有用信号。node
2、小波去噪的原理
Donoho提出的小波阀值去噪的基本思想是将信号经过小波变换（采用Mallat算法）后，信号产生的小波系数含有信号的重要信息，将信号经小波分解后小波系数较大，噪声的小波系数较小，而且噪声的小波系数要小于信号的小波系数，经过选取一个合适的阀值，大于阀值的小波系数被认为是有信号产生的，应予以保留，小于阀值的则认为是噪声产生的，置为零从而达到去噪的目的。python
从信号学的角度看 ,小波去噪是一个信号滤波的问题。尽管在很大程度上小波去噪能够当作是低通滤波 ,但因为在去噪后 ,还能成功地保留信号特征 ,因此在这一点上又优于传统的低通滤波器。因而可知 ,小波去噪其实是特征提取和低通滤波的综合 ,其流程图以下所示：算法
一个含噪的模型能够表示以下：app
其中 ,f( k)为有用信号,s(k)为含噪声信号,e(k)为噪声,ε为噪声系数的标准误差。less
假设，e(k)为高斯白噪声,一般状况下有用信号表现为低频部分或是一些比较平稳的信号,而噪声信号则表现为高频的信号,咱们对 s(k)信号进行小波分解的时候,则噪声部分一般包含在HL、LH、HH中,以下图所示，只要对HL、LH、HH做相应的小波系数处理,而后对信号进行重构便可以达到消噪的目的。函数
咱们能够看到，小波去噪的原理是比较简单类，相似以往咱们常见的低通滤波器的方法，可是因为小波去找保留了特征提取的部分，因此性能上是优于传统的去噪方法的。工具
3、小波去噪的基本方法
通常来讲， 一维信号的降噪过程能够分为 3个步骤性能
信号的小波分解。选择一个小波并肯定一个小波分解的层次N，而后对信号进行N层小波分解计算。spa
小波分解高频系数的阈值量化。对第1层到第N层的每一层高频系数（三个方向）， 选择一个阈值进行阈值量化处理．
这一步是最关键的一步，主要体如今阈值的选择与量化处理的过程，在每层阈值的选择上matlab提供了不少自适应的方法， 这里不一一介绍，量化处理方法主要有硬阈值量化与软阈值量化。下图是两者的区别：
上面左图是硬阈值量化，右图是软阈值量化。采用两种不一样的方法，达到的效果是，硬阈值方法能够很好地保留信号边缘等局部特征，软阈值处理相对要平滑，但会形成边缘模糊等失真现象。
信号的小波重构。根据小波分解的第 N层的低频系数和通过量化处理后的第1层到第N 层的高频系数，进行信号的小波重构。
小波阀值去噪的基本问题包括三个方面：小波基的选择，阀值的选择，阀值函数的选择。
（1）小波基的选择：一般咱们但愿所选取的小波知足如下条件：正交性、高消失矩、紧支性、对称性或反对称性。但事实上具备上述性质的小波是不可能存在的，由于小波是对称或反对称的只有Haar小波，而且高消失矩与紧支性是一对矛盾，因此在应用的时候通常选取具备紧支的小波以及根据信号的特征来选取较为合适的小波。
（2）阀值的选择：直接影响去噪效果的一个重要因素就是阀值的选取，不一样的阀值选取将有不一样的去噪效果。目前主要有通用阀值（VisuShrink）、SureShrink阀值、Minimax阀值、BayesShrink阀值等。
（3）阀值函数的选择：阀值函数是修正小波系数的规则，不一样的反之函数体现了不一样的处理小波系数的策略。最经常使用的阀值函数有两种：一种是硬阀值函数，另外一种是软阀值函数。还有一种介于软、硬阀值函数之间的Garrote函数。
另外，对于去噪效果好坏的评价，经常使用信号的信噪比（SNR）与估计信号同原始信号的均方根偏差（RMSE）来判断。
参考：
小波阀值去噪法基础 http://blog.sina.com.cn/s/blog_4d7c97a00101cib3.html
在python中使用小波分析进行阈值去噪声,使用pywt.threshold函数
#coding=gbk #使用小波分析进行阈值去噪声,使用pywt.threshold import pywt import numpy as np import pandas as pd import matplotlib.pyplot as plt import math data = np.linspace(1, 10, 10) print(data) # [ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] # pywt.threshold(data, value, mode, substitute) mode 模式有4种，soft, hard, greater, less; substitute是替换值 data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12) print(data_soft) # [12. 12. 12. 12. 12. 0. 1. 2. 3. 4.] 将小于6 的值设置为12， 大于等于6 的值所有减去6 data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12) print(data_hard) # [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12， 其他的值不变 data_greater = pywt.threshold(data, 6, 'greater', 12) print(data_greater) # [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12，大于等于阈值的值不变化 data_less = pywt.threshold(data, 6, 'less', 12) print(data_less) # [ 1. 2. 3. 4. 5. 6. 12. 12. 12. 12.] 将大于6 的值设置为12， 小于等于阈值的值不变
一、小波变换的理解
傅里叶变换——短时傅里叶变换——小波变换。
参考文献：如下两篇参考资料讲述得十分清楚，有助于理解小波变换。
但具体的数学角度阐述，请参考其余资料。
（1）知乎专栏：形象易懂讲解算法I——小波变换
https://zhuanlan.zhihu.com/p/22450818
（2）知乎专栏：傅里叶分析之掐死教程。
https://zhuanlan.zhihu.com/p/19763358
二、小波包分解
小波包是为了克服小波分解在高频段的频率分辨率较差，而在低频段的时间分辨率较差的问题的基础上而提出的。
它是一种更精细的信号分析的方法，提升了信号的时域分辨率。
下面是二者的对比图：
三、能量谱
基于小波包分解提取多尺度空间能量特征的原理是把不一样分解尺度上的信号能量求解出来，将这些能量值按尺度顺序排列成特征向量供识别使用。
20180510补充更新：具体计算公式以下所示，本文中未使用重构后的系数进行能量值计算，直接使用小波包分解后的系数，参考文献《基于小波包能量特征的滚动轴承故障监测方法 》。
四、Matlab代码
给出两部分代码，写成两个函数。一个是小波包分解与重构，另外一个是能量谱函数。
下载地址：https://download.csdn.net/download/ckzhb/10030651
代码名称：wavelet_packetdecomposition_reconstruct
function wpt= wavelet_packetdecomposition_reconstruct( x,n,wpname )
%% 对信号进行小波包分解，获得节点的小波包系数。而后对每一个节点系数进行重构。
% Decompose x at depth n with wpname wavelet packets.using Shannon entropy.
%
% x-input signal,列向量。
% n-the number of decomposition layers
% wpname-a particular wavelet.type:string.
%
%Author hubery_zhang
%Date 20170714
%%
wpt=wpdec(x,n,wpname);
% Plot wavelet packet tree (binary tree)
plot(wpt)
%% wavelet packet coefficients.default：use the front 4.
cfs0=wpcoef(wpt,[n 0]);
cfs1=wpcoef(wpt,[n 1]);
cfs2=wpcoef(wpt,[n 2]);
cfs3=wpcoef(wpt,[n 3]);
figure;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(cfs0);
title(['结点 ',num2str(n) ' 1',' 系数'])
subplot(5,1,3);
plot(cfs1);
title(['结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(cfs2);
title(['结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(cfs3);
title(['结点 ',num2str(n) ' 4',' 系数'])
%% reconstruct wavelet packet coefficients.
rex0=wprcoef(wpt,[n 0]);
rex1=wprcoef(wpt,[n 1]);
rex2=wprcoef(wpt,[n 2]);
rex3=wprcoef(wpt,[n 3]);
figure;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(rex0);
title(['重构结点 ',num2str(n) ' 1',' 系数'])
subplot(5,1,3);
plot(rex1);
title(['重构结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(rex2);
title(['重构结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(rex3);
title(['重构结点 ',num2str(n) ' 4',' 系数'])
end
代码名称：wavelet_energy_spectrum
function E = wavelet_energy_spectrum( wpt,n )
%% 计算每一层每个节点的能量
% wpt-wavelet packet tree
% n-第n层能量
%
% Author hubery_zhang
% Date 20170714
%%
% 求第n层第i个节点的系数
E(1:2^n )=0;
for i=1:2^n
E(i) = norm(wpcoef(wpt,[n,i-1]),2)^2; %20180604更新 原代码：E(i) = norm(wpcoef(wpt,[n,i-1]),2)
end
%求每一个节点的几率
E_total=sum(E);
for i=1:2^n
p_node(i)= 100*E(i)/E_total;
end
% E = wenergy(wpt); only get the last layer
figure;
x=1:2^n;
bar(x,p_node);
title(['第',num2str(n),'层']);
axis([0 2^n 0 100]);
xlabel('结点');
ylabel('能量百分比/%');
for j=1:2^n
text(x(j),p_node(i),num2str(p_node(j),'%0.2f'),...
'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
end
end
一维、二维离散卷积的计算方法：
1、 定义
离散信号f(n),g(n)的定义以下：
N-----为信号f(n)的长度
s(n)----为卷积结果序列,长度为len(f(n))+len(g(n))-1
以3个元素的信号为例：
f(n) = [1 2 3]; g(n) = [2 3 1];
s(0) = f(0)g(0-0) + f(1)g(0-1)+f(2)g(0-2) = 1*2 + 2*0 + 3*0 =2
s(1) = f(0)g(1-0) + f(1)g(1-1) + f(2)g(1-2) = 1*3 + 2*2 + 3*0 = 7
s(2) = f(0)g(2-0) + f(1)g(2-1) + f(2)g(2-2) =1*1 + 2*3 + 3*2=13
s(3) = f(0)g(3-0) + f(1)g(3-1) + f(2)g(3-2) =1*0 + 2*1 + 3*3=11
s(4) = f(0)g(4-0) + f(1)g(4-1) + f(2)g(4-2) =1*0 + 2*0 + 3*1=3
最终结果为:
s(n) = [2 7 13 11 3]
上述计算图示以下：
在数学里咱们知道f(-x)的图像是f(x)对y轴的反转
g(-m)就是把g(m)的序列反转，g(n-m)的意义是把g(-m)平移的n点：
如上图g(m)在信号处理中一般叫作滤波器或掩码，卷积至关于掩码g(m)反转后在信号f(n)上平移求和。Matlab计算卷积的函数为conv,
>> A = [1 2 3];
B = [2,3,1];
convD = conv(A,B)
convD =
2 7 13 11 3
相应的二维卷积定义以下：