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
相應的二維卷積定義以下：