小波包分解与重构
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
————————————————
版權聲明：本文為CSDN博主「hubery_zhang」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/ckzhb/article/details/78288847
?
?小波與小波包、小波包分解與信號重構、小波包能量特征提取? ?(Matlab 程序詳解）
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?-----暨 小波包分解后解決頻率大小分布重新排列問題
? ? ? ?本人當前對小波理解不是很深入，通過翻閱網絡他人博客，進行匯總總結，重新調試Matlab代碼，實現對小波與小波包、小波包分解與信號重構、小波包能量特征提取，供大家參考，后續將繼續更新！
? ? ?本人在分析信號的過程中發現，按照網上所述的小波包分解方法理解，獲取每層節點重構后信號頻率并不是按照（n,0）、(n,1)...順序依次由小到大排列的，經過進一步分析研究后發現，需要對節點進行重排序，具體操作見本文分析。
1.小波與小波包區別
? ? ? ? 工程應用中經常需要對一些非平穩信號進行，小波分析和小波包分析適合對非平穩信號分析，相比較小波分析，利用小波包分析可以對信號分析更加精細，小波包分析可以將時頻平面劃分的更為細致，對信號的高頻部分的分辨率要好于小波分析，可以根據信號的特征，自適應的選擇最佳小波基函數，比便更好的對信號進行分析，所以小波包分析應用更加廣泛。
? ? ? ? ? ? ? ? ? ? ? ? ?? ??
①小波分解
? ? ? ? ?小波變換只對信號的低頻部分做進一步分解，而對高頻部分也即信號的細節部分不再繼續分解，所以小波變換能夠很好地表征一大類以低頻信息為主要成分的信號，不能很好地分解和表示包含大量細節信息（細小邊緣或紋理）的信號，如非平穩機械振動信號、遙感圖象、地震信號和生物醫學信號等。
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
②小波包分解
? ? ? ?小波包變換既可以對低頻部分信號進行分解，也可以對高頻部分進行分解，而且這種分解既無冗余，也無疏漏，所以對包含大量中、高頻信息的信號能夠進行更好的時頻局部化分析。
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
2.小波包——小波包樹與時頻圖
小波包樹解讀：
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
?
? 以上即是小波包樹，其中節點的命名規則是從（1，0）開始，叫1號， （1，1）是2號………依此類推，（3，0）是7號，（3，7）是14號。 每個節點都有對應的小波包系數，這個系數決定了頻率的大小，也就是說頻率信息已經有了，但是時域信息在哪里呢？ 那就是 order。? 這個order就是這些節點的順序，也就是頻率的順序。
Matlab實例：
采樣頻率為1024Hz，采樣時間是1秒，有一個信號s是由頻率100和200Hz的正弦波混合的，我們用小波包來分解。
clear all
clc
fs=1024; %采樣頻率
f1=100; %信號的第一個頻率
f2=300; %信號第二個頻率
t=0:1/fs:1;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %生成混合信號
[tt]=wpdec(s,3,'dmey'); %小波包分解，3代表分解3層，'dmey'使用meyr小波
plot(tt) %畫小波包樹圖
wpviewcf(tt,1); %畫出時間頻率圖
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
主要解釋：
x軸，就是1024個點，對應1秒，每個點就代表1/1024秒。
y軸，顯示的數字對應于小波包樹中的節點，從下面開始，順序是 7號節點，8號，10號，9號，，，，11號節點，這個順序是這么排列的，這是小波包自動排列的。然后，y軸是頻率啊，怎么不是 100Hz和300Hz呢？我們的采樣頻率是1024Hz，根據采樣定理，奈奎斯特采樣頻率是512Hz，我們分解了3層，最后一層就是 2^3=8個頻率段，每個頻率段的頻率區間是 512/8=64Hz,看圖顏色重的地方一個是在8那里，一個在13那里，8是第二段，也就是 65-128Hz之間，13是第五段，也就是257-320Hz之間。這樣就說通了，正好這個原始信號只有兩個頻率段，一個100一個300。如果我們不是分解了3層，而是更多層，那么每個頻率段包含的頻率也就越窄，圖上有顏色的地方也會更細，也就是說更精細了，由于原始信號的頻率在整個1秒鐘內都沒有改變，所以有顏色的地方是一個橫線。（引用：http://www.cnblogs.com/welen/articles/5667217.html?）
3.小波包-----小波包分解系數
? ? ? ?在數值分析中，我們學過內積，內積的物理含義：兩個圖形的相似性，若兩個圖形完全正交，則內積為0，若兩個圖形完全一樣，則系數為1（相對值）。小波變換的實質是：原信號與小波基函數的相似性。小波系數就是小波基函數與原信號相似的系數。
? ? ? 連續小波變換：小波函數與原信號對應點相乘，再相加，得到對應點的小波變換系數，平移小波基函數，再計算小波函數與原信號對應點相乘，再相加，這樣就得到一系列的小波系數。對于離散小波變換（由于很多小波函數不是正交函數，因此需要一個尺度函數）所以，原信號函數可以分解成尺度函數和小波函數的線性組合，在這個函數中，尺度函數產生低頻部分，小波函數產生高頻部分。
4.小波包-----信號分解與重構（方法1）
該方法可以實現對任意節點系數選擇進行組合重構。
例1
有一個信號，變量名為wave，隨便找一個信號load進來就行了。
t=wpdec(wave,3,'dmey');
t2 = wpjoin(t,[3;4;5;6]);
sNod = read(t,'sizes',[3,4,5,6]);
cfs3 = zeros(sNod(1,:));
cfs4 = zeros(sNod(2,:));
cfs5 = zeros(sNod(3,:));
cfs6 = zeros(sNod(4,:));
t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);
wave2=wprec(t3);
第一行：將wave 用 dmey小波進行3層小波包分解，獲得一個小波包樹 t
第二行：將小波包樹的第二行的四個節點收起來，也就是讓第二行的節點變為樹的最底層節點。因為第一行中小波包樹的節點個數是第一層2個，第二層4個，第三層8個。現在t2就是將第三層的節點再聚合回第二層。
第三行：讀取第二層四個節點系數的size
第四~七行：將所有四個節點的小波包系數變為0
第八行：將四個節點的系數重組到t3小波樹中。
第九行：對t3小波樹進行重構，獲得信號wave2
? ? ? ?可以預見，因為我們把小波樹的節點系數都變為0了，所以信號也就全為0了。所以wave2是一個0向量。讀者可以自行plot一下wave和wave2看看。進一步，如果我們只聚合第二層中的某幾個節點，比如 4和5，即將第三行到第八行中節點 3 和節點 6的語句刪除或修改，那么意思就是將 4? 5節點的系數變為0，那么wave2肯定就不是0向量了。
例2
t=wpdec(wave,3,'dmey');
t2 = wpjoin(t,[3;4;5;6]);
cfs3=wpcoef(t,3);
cfs4=wpcoef(t,4);
cfs5=wpcoef(t,5);
cfs6=wpcoef(t,6);
t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);
wave2=wprec(t3);
解釋：
第一行：將wave 用 dmey小波進行3層小波包分解，獲得一個小波包樹 t
第二行：將小波包樹的第二行的四個節點收起來，也就是讓第二行的節點變為樹的最底層節點。
第三~六行：獲取四個節點的小波包系數 （小波包系數就是一個一維向量）
第七行：將四個節點的系數重組到t3小波樹中
第八行：對t3小波樹進行重構，獲得信號wave2
可以看出，該例子就是對一個小波包展開了，又原封不動的裝回去了，所以說 wave2和wave是一樣的。
注意，wpjoin命令在這里是必要的，因為write函數只能將最底層的節點寫進去。也就是說，如果我們將第三層的小波包系數進行修改的話，就不用wpjoin了，具體可以看例3
例3
t=wpdec(wave,3,'dmey');
cfs7=wpcoef(t,7);
cfs8=wpcoef(t,8);
cfs9=wpcoef(t,9);
cfs10=wpcoef(t,10);
cfs11=wpcoef(t,11);
cfs12=wpcoef(t,12);
cfs13=wpcoef(t,13);
cfs14=wpcoef(t,14);
t3=write(t,'cfs',7,cfs7,'cfs',8,cfs8,'cfs',9,cfs9,'cfs',10,cfs10,'cfs',11,cfs11,'cfs',...
12,cfs12,'cfs',13,cfs13,'cfs',14,cfs14);
y=wprec(t3);
該例子也是對一個小波包展開了，又原封不動的裝回去了，只不過這次是直接對第三層節點進行的。
5.小波包-----信號分解與重構（方法2）
該方法只能對某一節點信號系數分別進行重構，不能實現多個節點系數組合進行重構.
接下來進行舉例說明：
x_input=x_train(:,1,1); %輸入數據
plot(x_input);title('輸入信號時域圖像') %繪制輸入信號時域圖像
%% 查看頻譜范圍
x=x_input;
fs=128;
N=length(x); %采樣點個數
signalFFT=abs(fft(x,N));%真實的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
figure;plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');title('輸入信號的頻譜');grid on
接下來進行4層小波包分解
wpt=wpdec(x_input,3,'dmey'); %進行3層小波包分解
plot(wpt); %繪制小波包樹
接下來，我們查看第3層8個節點的頻譜分布（我的輸入信號采樣頻率是128，按照采樣定理小波包分解根節點（0，0）處的頻率應該為0-64Hz，按照這個計算（3，0）節點為0-8Hz，后面依次以8Hz為一個段遞增）
首先展示最初的重構方法（頻率排布順序混亂，常見理解錯誤）
for i=0:7
rex3(:,i+1)=wprcoef(wpt,[3 i]); %實現對節點小波節點進行重構
end
figure; %繪制第3層各個節點分別重構后信號的頻譜
for i=0:7
subplot(2,4,i+1);
x_sign=rex3(:,i+1);
N=length(x_sign); %采樣點個數
signalFFT=abs(fft(x_sign,N));%真實的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');grid on
axis([0 50 0 0.03]); title(['小波包第3層',num2str(i),'節點信號頻譜']);
end
繪制完后你會發現頻譜分布并不是按照之前理解的順序依次遞增排列，如下所示
那么問題出在哪里呢？經過仔細深入驗證后發現問題出在小波包節點的頻譜劃分“不是嚴格按照上述理解的順序排列的”（可能是一種格雷編碼或者其他），要解決這個問題我們需要對節點順序進行重新編碼排序。參考這里https://www.ilovematlab.cn/thread-122226-1-1.html?tdsourcetag=s_pctim_aiomsg
nodes=[7;8;9;10;11;12;13;14]; %第3層的節點號
ord=wpfrqord(nodes); ?%小波包系數重排，ord是重排后小波包系數索引構成的矩陣 如3層分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系數
注意：節點的命名規則是從（1，0）開始，叫1號， （1，1）是2號………依此類推，（3，0）是7號，（3，7）是14號。
那么我們來看看經過上面這段重排序運行后nodes_ord中順序是什么？
接下來我們再繪制一下第三層8個節點重構信號的頻譜如下
for i=1:8
rex3(:,i)=wprcoef(wpt,nodes_ord(i)); %實現對節點小波節點進行重構
end
figure; %繪制第3層各個節點分別重構后信號的頻譜
for i=0:7
subplot(2,4,i+1);
x_sign= rex3(:,i+1);
N=length(x_sign); %采樣點個數
signalFFT=abs(fft(x_sign,N));%真實的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');grid on
axis([0 50 0 0.03]); title(['小波包第3層',num2str(i),'節點信號頻譜']);
end
到這里為止不知道大家有沒有明白我要表達的意思，如果沒有明白可以再反復看理解，或者自己進行分解后觀察每個節點的頻譜后或許就理解了。
6.小波包分解------能量特征提取（方法1）
接下來繪制第3層各個頻段的能量占比
%% wavelet packet coefficients. 求取小波包分解的各個節點的小波包系數
cfs3_0=wpcoef(wpt,nodes_ord(1)); %對重排序后第3層0節點的小波包系數0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2)); %對重排序后第3層0節點的小波包系數8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3)); %對重排序后第3層0節點的小波包系數16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4)); %對重排序后第3層0節點的小波包系數24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5)); %對重排序后第3層0節點的小波包系數32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6)); %對重排序后第3層0節點的小波包系數40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7)); %對重排序后第3層0節點的小波包系數48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8)); %對重排序后第3層0節點的小波包系數56-64Hz
E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范數：就是norm(...,1)，即各元素絕對值之和；2-范數：就是norm(...,2)，即各元素平方和開根號；
E_cfs3_1=norm(cfs3_1,2)^2;
E_cfs3_2=norm(cfs3_2,2)^2;
E_cfs3_3=norm(cfs3_3,2)^2;
E_cfs3_4=norm(cfs3_4,2)^2;
E_cfs3_5=norm(cfs3_5,2)^2;
E_cfs3_6=norm(cfs3_6,2)^2;
E_cfs3_7=norm(cfs3_7,2)^2;
E_total=E_cfs3_0+E_cfs3_1+E_cfs3_2+E_cfs3_3+E_cfs3_4+E_cfs3_5+E_cfs3_6+E_cfs3_7;
p_node(0)= 100*E_cfs3_0/E_total; % 求得每個節點的占比
p_node(2)= 100*E_cfs3_1/E_total; % 求得每個節點的占比
p_node(3)= 100*E_cfs3_2/E_total; % 求得每個節點的占比
p_node(4)= 100*E_cfs3_3/E_total; % 求得每個節點的占比
p_node(5)= 100*E_cfs3_4/E_total; % 求得每個節點的占比
p_node(6)= 100*E_cfs3_5/E_total; % 求得每個節點的占比
p_node(7)= 100*E_cfs3_6/E_total; % 求得每個節點的占比
p_node(8)= 100*E_cfs3_7/E_total; % 求得每個節點的占比
figure;
x=1:8;
bar(x,p_node);
title('各個頻段能量所占的比例');
xlabel('頻率 Hz');
ylabel('能量百分比/%');
for j=1:8
text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...
'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
end
7.小波包分解------能量特征提取（方法2）
直接運行matlab自帶函數，如下
E = wenergy(wpt); ??%該函數只能對最后一層（底層）節點進行能量提取（引用：https://blog.csdn.net/it_beecoder/article/details/78668273）
這里對比一下，和方法1得到的圖形基本上是一致的，不同之處在于排列順序變了，方法2中的順序是按照小波包分解函數自動排列順序，方法1經過重排序后為按照頻率段遞增順序依次排列順序的
%% 繪制重構前各個頻段小波包系數
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信號');
subplot(4,1,2);
plot(cfs3_0);
title(['層數 ',num2str(3) ' 節點 0的小波0-8Hz',' 系數'])
subplot(4,1,3);
plot(cfs3_1);
title(['層數 ',num2str(3) ' 節點 1的小波8-16Hz',' 系數'])
subplot(4,1,4);
plot(cfs3_2);
title(['層數 ',num2str(3) ' 節點 2的小波16-24Hz',' 系數'])
%% 繪制重構后時域各個特征頻段的圖形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重構后0-8Hz頻段信號');
subplot(3,1,2);
plot(rex3(:,2));
title('重構后8-16Hz頻段信號')
subplot(3,1,3);
plot(rex3(:,3));
title('重構后16-24Hz頻段信號');
以上過程完整代碼
clear all;
x_input=x_train(:,1,1); %輸入數據
plot(x_input);title('輸入信號時域圖像') %繪制輸入信號時域圖像
x=x_input; % 查看頻譜范圍
fs=128;
N=length(x); %采樣點個數
signalFFT=abs(fft(x,N));%真實的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
figure;plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');title('輸入信號的頻譜');grid on
wpt=wpdec(x_input,3,'dmey'); %進行4層小波包分解
plot(wpt); %繪制小波包樹
%% 實現對節點順序按照頻率遞增進行重排序
% nodes=get(wpt,'tn'); %小波包分解系數 為什么ｎｏｄｅｓ是[7;8;9;10;11;12;13;14]
% N_cfs=length(nodes); %小波包系數個數
nodes=[7;8;9;10;11;12;13;14];
ord=wpfrqord(nodes); %小波包系數重排，ord是重排后小波包系數索引構成的矩陣 如3層分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系數
%% 實現對節點小波節點進行重構
for i=1:8
rex3(:,i)=wprcoef(wpt,nodes_ord(i));
end
%% 繪制第3層各個節點分別重構后信號的頻譜
figure;
for i=0:7
subplot(2,4,i+1);
x_sign= rex3(:,i+1);
N=length(x_sign); %采樣點個數
signalFFT=abs(fft(x_sign,N));%真實的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');grid on
axis([0 50 0 0.03]); title(['小波包第3層',num2str(i),'節點信號頻譜']);
end
%% wavelet packet coefficients. 求取小波包分解的各個頻段的小波包系數
cfs3_0=wpcoef(wpt,nodes_ord(1)); %對重排序后第3層0節點的小波包系數0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2)); %對重排序后第3層0節點的小波包系數8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3)); %對重排序后第3層0節點的小波包系數16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4)); %對重排序后第3層0節點的小波包系數24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5)); %對重排序后第3層0節點的小波包系數32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6)); %對重排序后第3層0節點的小波包系數40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7)); %對重排序后第3層0節點的小波包系數48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8)); %對重排序后第3層0節點的小波包系數56-64Hz
%% 求取小波包分解的各個頻段的能量
E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范數：就是norm(...,1)，即各元素絕對值之和；2-范數：就是norm(...,2)，即各元素平方和開根號；
E_cfs3_1=norm(cfs3_1,2)^2;
E_cfs3_2=norm(cfs3_2,2)^2;
E_cfs3_3=norm(cfs3_3,2)^2;
E_cfs3_4=norm(cfs3_4,2)^2;
E_cfs3_5=norm(cfs3_5,2)^2;
E_cfs3_6=norm(cfs3_6,2)^2;
E_cfs3_7=norm(cfs3_7,2)^2;
E_total=E_cfs3_0+E_cfs3_1+E_cfs3_2+E_cfs3_3+E_cfs3_4+E_cfs3_5+E_cfs3_6+E_cfs3_7;
p_node(0)= 100*E_cfs3_0/E_total; % 求得每個節點的占比
p_node(2)= 100*E_cfs3_1/E_total; % 求得每個節點的占比
p_node(3)= 100*E_cfs3_2/E_total; % 求得每個節點的占比
p_node(4)= 100*E_cfs3_3/E_total; % 求得每個節點的占比
p_node(5)= 100*E_cfs3_4/E_total; % 求得每個節點的占比
p_node(6)= 100*E_cfs3_5/E_total; % 求得每個節點的占比
p_node(7)= 100*E_cfs3_6/E_total; % 求得每個節點的占比
p_node(8)= 100*E_cfs3_7/E_total; % 求得每個節點的占比
figure;
x=1:8;
bar(x,p_node);
title('各個頻段能量所占的比例');
xlabel('頻率 Hz');
ylabel('能量百分比/%');
for j=1:8
text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...
'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
end
% E = wenergy(wpt); %求取各個節點能量
%% 繪制重構前各個特征頻段小波包系數的圖形
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信號');
subplot(4,1,2);
plot(cfs3_0);
title(['層數 ',num2str(3) ' 節點 0的小波0-8Hz',' 系數'])
subplot(4,1,3);
plot(cfs3_1);
title(['層數 ',num2str(3) ' 節點 1的小波8-16Hz',' 系數'])
subplot(4,1,4);
plot(cfs3_2);
title(['層數 ',num2str(3) ' 節點 2的小波16-24Hz',' 系數'])
%% 繪制重構后時域各個特征頻段的圖形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重構后0-8Hz頻段信號');
subplot(3,1,2);
plot(rex3(:,2));
title('重構后8-16Hz頻段信號')
subplot(3,1,3);
plot(rex3(:,3));
title('重構后16-24Hz頻段信號');
8.小波----常見基函數
? ? ?與標準的傅里葉變換相比，小波分析中使用到的小波函數具有不唯一性，即小波函數 具有多樣性。小波分析在工程應用中，一個十分重要的問題就是最優小波基的選擇問題，因為用不同的小波基分析同一個問題會產生不同的結果。目前我們主要是通過用小波分析方法處理信號的結果與理論結果的誤差來判定小波基的好壞，由此決定小波基。
???? 常用小波基有Haar小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet小波、Meyer小波等。
?????1.Haar小波Haar函數是小波分析中最早用到的一個具有緊支撐的正交小波函數，也是最簡單的一個小波函數，它是支撐域在范圍內的單個矩形波。Haar函數的定義如下：Haar小波在時域上是不連續的，所以作為基本小波性能不是特別好。但它也有自己的優點:計算簡單。不但與正交，而且與自己的整數位移正交，因此，在的多分辨率系統中，Haar小波構成一組最簡單的正交歸一的小波族。
[phi,g1,xval] = wavefun('haar',20);
subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t')；
title('haar 時域');
g2=fft(g1);
g3=abs(g2);
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f')title('haar 頻域')；
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
2.Daubechies(dbN)小波Daubechies小波是世界著名的小波分析學者Inrid·Daubechies構造的小波函數，簡寫為dbN，N是小波的階數。小波和尺度函數中的支撐區為的消失矩為。除（Harr小波）外，dbN不具有對稱性（即非線性相位）。除（Harr小波）外，dbN沒有明確的表達式，但轉換函數h的平方模是明確的:令，其中為二項式的系數，則有其中：Daubechies小波具有以下特點：在時域是有限支撐的，即長度有限。在頻域在處有N階零點。和它的整數位移正交歸一，即。小波函數可以由所謂“尺度函數”求出來。尺度函數為低通函數，長度有限，支撐域在的范圍內。
db4的時域和頻域波形：
[phi,g1,xval] = wavefun('db4',10);
subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t')title('db4 時域');
g2=fft(g1);
g3=abs(g2);
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f')title('db4 頻域')；
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
Daubechies小波常用來分解和重構信號，作為濾波器使用:
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db4'); %計算該小波的4個濾波器
subplot(2,2,1);
stem(Lo_D,'LineWidth',2);
title('分解低通濾波器');
subplot(2,2,2);
stem(Hi_D,'LineWidth',2);
title('分解高通濾波器');
subplot(2,2,3);
stem(Lo_R,'LineWidth',2);
title('重構低通濾波器');
subplot(2,2,4);
stem(Hi_R,'LineWidth',2);
title('重構高通濾波器');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
? 3.Mexican Hat(mexh)小波Mexican Hat函數為Gauss函數的二階導數：因為它的形狀像墨西哥帽的截面，所以也稱為墨西哥帽函數。Mexihat小波的時域和頻域波形：
Mexihat小波的時域和頻域波形：
d=-6; h=6; n=100;
[g1,x]=mexihat(d,h,n);
subplot(2,1,1);
plot(x,g1,'LineWidth',2);
xlabel('t');
title('Mexihat 時域');
g2=fft(g1);
g3=(abs(g2));
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 頻域');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
Mexihat小波的特點：在時間域與頻率域都有很好的局部化，并且滿足。不存在尺度函數，所以Mexihat小波函數不具有正交性。
?4.Morlet小波它是高斯包絡下的單頻率副正弦函數：其中C是重構時的歸一化常數。Morlet小波沒有尺度函數，而且是非正交分解。Morlet小波的時域和頻域波形圖：
d=-6; h=6; n=100;
[g1,x]=morlet(d,h,n);
subplot(2,1,1);
plot(x,g1,'LineWidth',2);
xlabel('t');
title('morlet 時域');
g2=fft(g1);
g3=(abs(g2));
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 頻域');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
注:
獲取小波系數的兩個函數理解：
Wpcoef?是求解某個節點的小波包系數，數據長度是1/(2^n)（分解n層的話），其實就是將原始信號化成2^N段，每段的長度是相等的且比原信號短
wprcoef是把某個節點的小波包系數重構，得到的是和原信號一樣長度的信號。
傅里葉變換與小波變換理解參考:
https://zhuanlan.zhihu.com/p/22450818
https://zhuanlan.zhihu.com/p/19763358
https://www.jianshu.com/p/5b5c160c7e3a
總結
- 上一篇： Vscode下中文乱码问题
- 下一篇： PS剪裁工具详解