一、前言
在現實生活和工作中，噪聲無處不在，在許多領域中，如天文、醫學圖像和計算機視覺方面收集到的數據常常是含有噪聲的。噪聲可能來自獲取數據的過程，也可能來自環境影響。由於種種原因,總會存在噪聲，噪聲的存在往往會掩蓋信號本身所要表現的信息，所以在實際的信號處理中，常常需要對信號進行預處理，而預處理最主要的一個步驟就是降噪。
小波分析是近年來發展起來的一種新的信號處理工具，這種方法源於傅立葉分析，小波（wavelet），即小區域的波，僅僅在非常有限的一段區間有非零值，而不是像正弦波和余弦波那樣無始無終。小波可以沿時間軸前后平移，也可按比例伸展和壓縮以獲取低頻和高頻小波，構造好的小波函數可以用於濾波或壓縮信號，從而可以提取出已含噪聲信號中的有用信號。
二、小波去噪的原理
Donoho提出的小波閥值去噪的基本思想是將信號通過小波變換（采用Mallat算法）后，信號產生的小波系數含有信號的重要信息，將信號經小波分解后小波系數較大，噪聲的小波系數較小，並且噪聲的小波系數要小於信號的小波系數，通過選取一個合適的閥值，大於閥值的小波系數被認為是有信號產生的，應予以保留，小於閥值的則認為是噪聲產生的，置為零從而達到去噪的目的。
從信號學的角度看 ,小波去噪是一個信號濾波的問題。盡管在很大程度上小波去噪可以看成是低通濾波 ,但由於在去噪后 ,還能成功地保留信號特征 ,所以在這一點上又優於傳統的低通濾波器。由此可見 ,小波去噪實際上是特征提取和低通濾波的綜合 ,其流程圖如下所示：
一個含噪的模型可以表示如下：
其中 ,f( k)為有用信號,s(k)為含噪聲信號,e(k)為噪聲,ε為噪聲系數的標准偏差。
假設，e(k)為高斯白噪聲,通常情況下有用信號表現為低頻部分或是一些比較平穩的信號,而噪聲信號則表現為高頻的信號,我們對 s(k)信號進行小波分解的時候,則噪聲部分通常包含在HL、LH、HH中,如下圖所示，只要對HL、LH、HH作相應的小波系數處理,然后對信號進行重構即可以達到消噪的目的。
我們可以看到，小波去噪的原理是比較簡單類，類似以往我們常見的低通濾波器的方法，但是由於小波去找保留了特征提取的部分，所以性能上是優於傳統的去噪方法的。
三、小波去噪的基本方法
一般來說， 一維信號的降噪過程可以分為 3個步驟
信號的小波分解。選擇一個小波並確定一個小波分解的層次N，然后對信號進行N層小波分解計算。
小波分解高頻系數的閾值量化。對第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、小波變換的理解
傅里葉變換——短時傅里葉變換——小波變換。
參考文獻：以下兩篇參考資料講述得十分清楚，有助於理解小波變換。
但具體的數學角度闡述，請參考其他資料。
（1）知乎專欄：形象易懂講解算法I——小波變換
https://zhuanlan.zhihu.com/p/22450818
（2）知乎專欄：傅里葉分析之掐死教程。
https://zhuanlan.zhihu.com/p/19763358
2、小波包分解
小波包是為了克服小波分解在高頻段的頻率分辨率較差，而在低頻段的時間分辨率較差的問題的基礎上而提出的。
它是一種更精細的信號分析的方法，提高了信號的時域分辨率。
下面是兩者的對比圖：
3、能量譜
基於小波包分解提取多尺度空間能量特征的原理是把不同分解尺度上的信號能量求解出來，將這些能量值按尺度順序排列成特征向量供識別使用。
20180510補充更新：具體計算公式如下所示，本文中未使用重構后的系數進行能量值計算，直接使用小波包分解后的系數，參考文獻《基於小波包能量特征的滾動軸承故障監測方法 》。
4、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一維、二維離散卷積的計算方法：
一、 定義
離散信號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
相應的二維卷積定義如下：
