小波与小波包、小波包分解与信号重构、小波包能量特征提取 暨 小波包分解后实现按频率大小分布重新排列(Matlab 程序详解）
小波與小波包、小波包分解與信號(hào)重構(gòu)、小波包能量特征提取? ?(Matlab 程序詳解）
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?-----暨 小波包分解后解決頻率大小分布重新排列問題
? ? ? ?本人當(dāng)前對(duì)小波理解不是很深入，通過翻閱網(wǎng)絡(luò)他人博客，進(jìn)行匯總總結(jié)，重新調(diào)試Matlab代碼，實(shí)現(xiàn)對(duì)小波與小波包、小波包分解與信號(hào)重構(gòu)、小波包能量特征提取，供大家參考，后續(xù)將繼續(xù)更新！
? ? ?本人在分析信號(hào)的過程中發(fā)現(xiàn)，按照網(wǎng)上所述的小波包分解方法理解，獲取每層節(jié)點(diǎn)重構(gòu)后信號(hào)頻率并不是按照（n,0）、(n,1)...順序依次由小到大排列的，經(jīng)過進(jìn)一步分析研究后發(fā)現(xiàn)，需要對(duì)節(jié)點(diǎn)進(jìn)行重排序，具體操作見本文分析。
1.小波與小波包區(qū)別
? ? ? ? 工程應(yīng)用中經(jīng)常需要對(duì)一些非平穩(wěn)信號(hào)進(jìn)行，小波分析和小波包分析適合對(duì)非平穩(wěn)信號(hào)分析，相比較小波分析，利用小波包分析可以對(duì)信號(hào)分析更加精細(xì)，小波包分析可以將時(shí)頻平面劃分的更為細(xì)致，對(duì)信號(hào)的高頻部分的分辨率要好于小波分析，可以根據(jù)信號(hào)的特征，自適應(yīng)的選擇最佳小波基函數(shù)，比便更好的對(duì)信號(hào)進(jìn)行分析，所以小波包分析應(yīng)用更加廣泛。
? ? ? ? ? ? ? ? ? ? ? ? ?? ??
①小波分解
? ? ? ? ?小波變換只對(duì)信號(hào)的低頻部分做進(jìn)一步分解，而對(duì)高頻部分也即信號(hào)的細(xì)節(jié)部分不再繼續(xù)分解，所以小波變換能夠很好地表征一大類以低頻信息為主要成分的信號(hào)，不能很好地分解和表示包含大量細(xì)節(jié)信息（細(xì)小邊緣或紋理）的信號(hào)，如非平穩(wěn)機(jī)械振動(dòng)信號(hào)、遙感圖象、地震信號(hào)和生物醫(yī)學(xué)信號(hào)等。
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
②小波包分解
? ? ? ?小波包變換既可以對(duì)低頻部分信號(hào)進(jìn)行分解，也可以對(duì)高頻部分進(jìn)行分解，而且這種分解既無冗余，也無疏漏，所以對(duì)包含大量中、高頻信息的信號(hào)能夠進(jìn)行更好的時(shí)頻局部化分析。
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
2.小波包——小波包樹與時(shí)頻圖
小波包樹解讀：
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
?
? 以上即是小波包樹，其中節(jié)點(diǎn)的命名規(guī)則是從（1，0）開始，叫1號(hào)， （1，1）是2號(hào)………依此類推，（3，0）是7號(hào)，（3，7）是14號(hào)。 每個(gè)節(jié)點(diǎn)都有對(duì)應(yīng)的小波包系數(shù)，這個(gè)系數(shù)決定了頻率的大小，也就是說頻率信息已經(jīng)有了，但是時(shí)域信息在哪里呢？ 那就是 order。? 這個(gè)order就是這些節(jié)點(diǎn)的順序，也就是頻率的順序。
Matlab實(shí)例：
采樣頻率為1024Hz，采樣時(shí)間是1秒，有一個(gè)信號(hào)s是由頻率100和200Hz的正弦波混合的，我們用小波包來分解。
?clear all
clc
fs=1024; %采樣頻率
f1=100; %信號(hào)的第一個(gè)頻率
f2=300; %信號(hào)第二個(gè)頻率
t=0:1/fs:1;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %生成混合信號(hào)
[tt]=wpdec(s,3,'dmey'); %小波包分解，3代表分解3層，'dmey'使用meyr小波
plot(tt) %畫小波包樹圖
wpviewcf(tt,1); %畫出時(shí)間頻率圖
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
主要解釋：
x軸，就是1024個(gè)點(diǎn)，對(duì)應(yīng)1秒，每個(gè)點(diǎn)就代表1/1024秒。
y軸，顯示的數(shù)字對(duì)應(yīng)于小波包樹中的節(jié)點(diǎn)，從下面開始，順序是 7號(hào)節(jié)點(diǎn)，8號(hào)，10號(hào)，9號(hào)，，，，11號(hào)節(jié)點(diǎn)，這個(gè)順序是這么排列的，這是小波包自動(dòng)排列的。然后，y軸是頻率啊，怎么不是 100Hz和300Hz呢？我們的采樣頻率是1024Hz，根據(jù)采樣定理，奈奎斯特采樣頻率是512Hz，我們分解了3層，最后一層就是 2^3=8個(gè)頻率段，每個(gè)頻率段的頻率區(qū)間是 512/8=64Hz,看圖顏色重的地方一個(gè)是在8那里，一個(gè)在13那里，8是第二段，也就是 65-128Hz之間，13是第五段，也就是257-320Hz之間。這樣就說通了，正好這個(gè)原始信號(hào)只有兩個(gè)頻率段，一個(gè)100一個(gè)300。如果我們不是分解了3層，而是更多層，那么每個(gè)頻率段包含的頻率也就越窄，圖上有顏色的地方也會(huì)更細(xì)，也就是說更精細(xì)了，由于原始信號(hào)的頻率在整個(gè)1秒鐘內(nèi)都沒有改變，所以有顏色的地方是一個(gè)橫線。（引用：http://www.cnblogs.com/welen/articles/5667217.html?）
3.小波包-----小波包分解系數(shù)
? ? ? ?在數(shù)值分析中，我們學(xué)過內(nèi)積，內(nèi)積的物理含義：兩個(gè)圖形的相似性，若兩個(gè)圖形完全正交，則內(nèi)積為0，若兩個(gè)圖形完全一樣，則系數(shù)為1（相對(duì)值）。小波變換的實(shí)質(zhì)是：原信號(hào)與小波基函數(shù)的相似性。小波系數(shù)就是小波基函數(shù)與原信號(hào)相似的系數(shù)。
? ? ? 連續(xù)小波變換：小波函數(shù)與原信號(hào)對(duì)應(yīng)點(diǎn)相乘，再相加，得到對(duì)應(yīng)點(diǎn)的小波變換系數(shù)，平移小波基函數(shù)，再計(jì)算小波函數(shù)與原信號(hào)對(duì)應(yīng)點(diǎn)相乘，再相加，這樣就得到一系列的小波系數(shù)。對(duì)于離散小波變換（由于很多小波函數(shù)不是正交函數(shù)，因此需要一個(gè)尺度函數(shù)）所以，原信號(hào)函數(shù)可以分解成尺度函數(shù)和小波函數(shù)的線性組合，在這個(gè)函數(shù)中，尺度函數(shù)產(chǎn)生低頻部分，小波函數(shù)產(chǎn)生高頻部分。
4.小波包-----信號(hào)分解與重構(gòu)（方法1）
該方法可以實(shí)現(xiàn)對(duì)任意節(jié)點(diǎn)系數(shù)選擇進(jìn)行組合重構(gòu)。
例1
有一個(gè)信號(hào)，變量名為wave，隨便找一個(gè)信號(hào)load進(jìn)來就行了。
?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小波進(jìn)行3層小波包分解，獲得一個(gè)小波包樹 t
第二行：將小波包樹的第二行的四個(gè)節(jié)點(diǎn)收起來，也就是讓第二行的節(jié)點(diǎn)變?yōu)闃涞淖畹讓庸?jié)點(diǎn)。因?yàn)榈谝恍兄行〔ò鼧涞墓?jié)點(diǎn)個(gè)數(shù)是第一層2個(gè)，第二層4個(gè)，第三層8個(gè)。現(xiàn)在t2就是將第三層的節(jié)點(diǎn)再聚合回第二層。
第三行：讀取第二層四個(gè)節(jié)點(diǎn)系數(shù)的size
第四~七行：將所有四個(gè)節(jié)點(diǎn)的小波包系數(shù)變?yōu)?
第八行：將四個(gè)節(jié)點(diǎn)的系數(shù)重組到t3小波樹中。
第九行：對(duì)t3小波樹進(jìn)行重構(gòu)，獲得信號(hào)wave2
? ? ? ?可以預(yù)見，因?yàn)槲覀儼研〔涞墓?jié)點(diǎn)系數(shù)都變?yōu)?了，所以信號(hào)也就全為0了。所以wave2是一個(gè)0向量。讀者可以自行plot一下wave和wave2看看。進(jìn)一步，如果我們只聚合第二層中的某幾個(gè)節(jié)點(diǎn)，比如 4和5，即將第三行到第八行中節(jié)點(diǎn) 3 和節(jié)點(diǎn) 6的語句刪除或修改，那么意思就是將 4? 5節(jié)點(diǎn)的系數(shù)變?yōu)?，那么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小波進(jìn)行3層小波包分解，獲得一個(gè)小波包樹 t
第二行：將小波包樹的第二行的四個(gè)節(jié)點(diǎn)收起來，也就是讓第二行的節(jié)點(diǎn)變?yōu)闃涞淖畹讓庸?jié)點(diǎn)。
第三~六行：獲取四個(gè)節(jié)點(diǎn)的小波包系數(shù) （小波包系數(shù)就是一個(gè)一維向量）
第七行：將四個(gè)節(jié)點(diǎn)的系數(shù)重組到t3小波樹中
第八行：對(duì)t3小波樹進(jìn)行重構(gòu)，獲得信號(hào)wave2
可以看出，該例子就是對(duì)一個(gè)小波包展開了，又原封不動(dòng)的裝回去了，所以說 wave2和wave是一樣的。
注意，wpjoin命令在這里是必要的，因?yàn)閣rite函數(shù)只能將最底層的節(jié)點(diǎn)寫進(jìn)去。也就是說，如果我們將第三層的小波包系數(shù)進(jìn)行修改的話，就不用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);
該例子也是對(duì)一個(gè)小波包展開了，又原封不動(dòng)的裝回去了，只不過這次是直接對(duì)第三層節(jié)點(diǎn)進(jìn)行的。
5.小波包-----信號(hào)分解與重構(gòu)（方法2）
該方法只能對(duì)某一節(jié)點(diǎn)信號(hào)系數(shù)分別進(jìn)行重構(gòu)，不能實(shí)現(xiàn)多個(gè)節(jié)點(diǎn)系數(shù)組合進(jìn)行重構(gòu).
接下來進(jìn)行舉例說明：
?x_input=x_train(:,1,1); %輸入數(shù)據(jù)
plot(x_input);title('輸入信號(hào)時(shí)域圖像') %繪制輸入信號(hào)時(shí)域圖像
%% 查看頻譜范圍
x=x_input;
fs=128;
N=length(x); %采樣點(diǎn)個(gè)數(shù)
signalFFT=abs(fft(x,N));%真實(shí)的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
figure;plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');title('輸入信號(hào)的頻譜');grid on
接下來進(jìn)行4層小波包分解
?wpt=wpdec(x_input,3,'dmey'); %進(jìn)行3層小波包分解
plot(wpt); %繪制小波包樹
接下來，我們查看第3層8個(gè)節(jié)點(diǎn)的頻譜分布（我的輸入信號(hào)采樣頻率是128，按照采樣定理小波包分解根節(jié)點(diǎn)（0，0）處的頻率應(yīng)該為0-64Hz，按照這個(gè)計(jì)算（3，0）節(jié)點(diǎn)為0-8Hz，后面依次以8Hz為一個(gè)段遞增）
首先展示最初的重構(gòu)方法（頻率排布順序混亂，常見理解錯(cuò)誤）
?for i=0:7
rex3(:,i+1)=wprcoef(wpt,[3 i]); %實(shí)現(xiàn)對(duì)節(jié)點(diǎn)小波節(jié)點(diǎn)進(jìn)行重構(gòu)
end
figure; %繪制第3層各個(gè)節(jié)點(diǎn)分別重構(gòu)后信號(hào)的頻譜
for i=0:7
subplot(2,4,i+1);
x_sign=rex3(:,i+1);
N=length(x_sign); %采樣點(diǎn)個(gè)數(shù)
signalFFT=abs(fft(x_sign,N));%真實(shí)的幅值
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),'節(jié)點(diǎn)信號(hào)頻譜']);
end
繪制完后你會(huì)發(fā)現(xiàn)頻譜分布并不是按照之前理解的順序依次遞增排列，如下所示
那么問題出在哪里呢？經(jīng)過仔細(xì)深入驗(yàn)證后發(fā)現(xiàn)問題出在小波包節(jié)點(diǎn)的頻譜劃分“不是嚴(yán)格按照上述理解的順序排列的”（可能是一種格雷編碼或者其他），要解決這個(gè)問題我們需要對(duì)節(jié)點(diǎn)順序進(jìn)行重新編碼排序。參考這里https://www.ilovematlab.cn/thread-122226-1-1.html?tdsourcetag=s_pctim_aiomsg
?nodes=[7;8;9;10;11;12;13;14]; %第3層的節(jié)點(diǎn)號(hào)
ord=wpfrqord(nodes); ?%小波包系數(shù)重排，ord是重排后小波包系數(shù)索引構(gòu)成的矩陣 如3層分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系數(shù)
注意：節(jié)點(diǎn)的命名規(guī)則是從（1，0）開始，叫1號(hào)， （1，1）是2號(hào)………依此類推，（3，0）是7號(hào)，（3，7）是14號(hào)。
那么我們來看看經(jīng)過上面這段重排序運(yùn)行后nodes_ord中順序是什么？
接下來我們?cè)倮L制一下第三層8個(gè)節(jié)點(diǎn)重構(gòu)信號(hào)的頻譜如下
?for i=1:8
rex3(:,i)=wprcoef(wpt,nodes_ord(i)); %實(shí)現(xiàn)對(duì)節(jié)點(diǎn)小波節(jié)點(diǎn)進(jìn)行重構(gòu)
end
figure; %繪制第3層各個(gè)節(jié)點(diǎn)分別重構(gòu)后信號(hào)的頻譜
for i=0:7
subplot(2,4,i+1);
x_sign= rex3(:,i+1);
N=length(x_sign); %采樣點(diǎn)個(gè)數(shù)
signalFFT=abs(fft(x_sign,N));%真實(shí)的幅值
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),'節(jié)點(diǎn)信號(hào)頻譜']);
end
到這里為止不知道大家有沒有明白我要表達(dá)的意思，如果沒有明白可以再反復(fù)看理解，或者自己進(jìn)行分解后觀察每個(gè)節(jié)點(diǎn)的頻譜后或許就理解了。
6.小波包分解------能量特征提取（方法1）
接下來繪制第3層各個(gè)頻段的能量占比
?%% wavelet packet coefficients. 求取小波包分解的各個(gè)節(jié)點(diǎn)的小波包系數(shù)
cfs3_0=wpcoef(wpt,nodes_ord(1)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)56-64Hz
E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范數(shù)：就是norm(...,1)，即各元素絕對(duì)值之和；2-范數(shù)：就是norm(...,2)，即各元素平方和開根號(hào)；
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; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(2)= 100*E_cfs3_1/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(3)= 100*E_cfs3_2/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(4)= 100*E_cfs3_3/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(5)= 100*E_cfs3_4/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(6)= 100*E_cfs3_5/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(7)= 100*E_cfs3_6/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(8)= 100*E_cfs3_7/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
figure;
x=1:8;
bar(x,p_node);
title('各個(gè)頻段能量所占的比例');
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）
直接運(yùn)行matlab自帶函數(shù)，如下
E = wenergy(wpt); ??%該函數(shù)只能對(duì)最后一層（底層）節(jié)點(diǎn)進(jìn)行能量提取（引用：https://blog.csdn.net/it_beecoder/article/details/78668273）
這里對(duì)比一下，和方法1得到的圖形基本上是一致的，不同之處在于排列順序變了，方法2中的順序是按照小波包分解函數(shù)自動(dòng)排列順序，方法1經(jīng)過重排序后為按照頻率段遞增順序依次排列順序的
?%% 繪制重構(gòu)前各個(gè)頻段小波包系數(shù)
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信號(hào)');
subplot(4,1,2);
plot(cfs3_0);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 0的小波0-8Hz',' 系數(shù)'])
subplot(4,1,3);
plot(cfs3_1);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 1的小波8-16Hz',' 系數(shù)'])
subplot(4,1,4);
plot(cfs3_2);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 2的小波16-24Hz',' 系數(shù)'])
%% 繪制重構(gòu)后時(shí)域各個(gè)特征頻段的圖形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重構(gòu)后0-8Hz頻段信號(hào)');
subplot(3,1,2);
plot(rex3(:,2));
title('重構(gòu)后8-16Hz頻段信號(hào)')
subplot(3,1,3);
plot(rex3(:,3));
title('重構(gòu)后16-24Hz頻段信號(hào)');
以上過程完整代碼
?clear all;
x_input=x_train(:,1,1); %輸入數(shù)據(jù)
plot(x_input);title('輸入信號(hào)時(shí)域圖像') %繪制輸入信號(hào)時(shí)域圖像
x=x_input; % 查看頻譜范圍
fs=128;
N=length(x); %采樣點(diǎn)個(gè)數(shù)
signalFFT=abs(fft(x,N));%真實(shí)的幅值
Y=2*signalFFT/N;
f=(0:N/2)*(fs/N);
figure;plot(f,Y(1:N/2+1));
ylabel('amp'); xlabel('frequency');title('輸入信號(hào)的頻譜');grid on
wpt=wpdec(x_input,3,'dmey'); %進(jìn)行4層小波包分解
plot(wpt); %繪制小波包樹
%% 實(shí)現(xiàn)對(duì)節(jié)點(diǎn)順序按照頻率遞增進(jìn)行重排序
% nodes=get(wpt,'tn'); %小波包分解系數(shù) 為什么ｎｏｄｅｓ是[7;8;9;10;11;12;13;14]
% N_cfs=length(nodes); %小波包系數(shù)個(gè)數(shù)
nodes=[7;8;9;10;11;12;13;14];
ord=wpfrqord(nodes); %小波包系數(shù)重排，ord是重排后小波包系數(shù)索引構(gòu)成的矩陣 如3層分解的[1;2;4;3;7;8;6;5]
nodes_ord=nodes(ord); %重排后的小波系數(shù)
%% 實(shí)現(xiàn)對(duì)節(jié)點(diǎn)小波節(jié)點(diǎn)進(jìn)行重構(gòu)
for i=1:8
rex3(:,i)=wprcoef(wpt,nodes_ord(i));
end
%% 繪制第3層各個(gè)節(jié)點(diǎn)分別重構(gòu)后信號(hào)的頻譜
figure;
for i=0:7
subplot(2,4,i+1);
x_sign= rex3(:,i+1);
N=length(x_sign); %采樣點(diǎn)個(gè)數(shù)
signalFFT=abs(fft(x_sign,N));%真實(shí)的幅值
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),'節(jié)點(diǎn)信號(hào)頻譜']);
end
%% wavelet packet coefficients. 求取小波包分解的各個(gè)頻段的小波包系數(shù)
cfs3_0=wpcoef(wpt,nodes_ord(1)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)0-8Hz
cfs3_1=wpcoef(wpt,nodes_ord(2)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)8-16Hz
cfs3_2=wpcoef(wpt,nodes_ord(3)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)16-24Hz
cfs3_3=wpcoef(wpt,nodes_ord(4)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)24-32Hz
cfs3_4=wpcoef(wpt,nodes_ord(5)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)32-40Hz
cfs3_5=wpcoef(wpt,nodes_ord(6)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)40-48Hz
cfs3_6=wpcoef(wpt,nodes_ord(7)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)48-56Hz
cfs3_7=wpcoef(wpt,nodes_ord(8)); %對(duì)重排序后第3層0節(jié)點(diǎn)的小波包系數(shù)56-64Hz
%% 求取小波包分解的各個(gè)頻段的能量
E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范數(shù)：就是norm(...,1)，即各元素絕對(duì)值之和；2-范數(shù)：就是norm(...,2)，即各元素平方和開根號(hào)；
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; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(2)= 100*E_cfs3_1/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(3)= 100*E_cfs3_2/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(4)= 100*E_cfs3_3/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(5)= 100*E_cfs3_4/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(6)= 100*E_cfs3_5/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(7)= 100*E_cfs3_6/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
p_node(8)= 100*E_cfs3_7/E_total; % 求得每個(gè)節(jié)點(diǎn)的占比
figure;
x=1:8;
bar(x,p_node);
title('各個(gè)頻段能量所占的比例');
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); %求取各個(gè)節(jié)點(diǎn)能量
%% 繪制重構(gòu)前各個(gè)特征頻段小波包系數(shù)的圖形
figure(1);
subplot(4,1,1);
plot(x_input);
title('原始信號(hào)');
subplot(4,1,2);
plot(cfs3_0);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 0的小波0-8Hz',' 系數(shù)'])
subplot(4,1,3);
plot(cfs3_1);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 1的小波8-16Hz',' 系數(shù)'])
subplot(4,1,4);
plot(cfs3_2);
title(['層數(shù) ',num2str(3) ' 節(jié)點(diǎn) 2的小波16-24Hz',' 系數(shù)'])
%% 繪制重構(gòu)后時(shí)域各個(gè)特征頻段的圖形
figure(3);
subplot(3,1,1);
plot(rex3(:,1));
title('重構(gòu)后0-8Hz頻段信號(hào)');
subplot(3,1,2);
plot(rex3(:,2));
title('重構(gòu)后8-16Hz頻段信號(hào)')
subplot(3,1,3);
plot(rex3(:,3));
title('重構(gòu)后16-24Hz頻段信號(hào)');
8.小波----常見基函數(shù)
? ? ?與標(biāo)準(zhǔn)的傅里葉變換相比，小波分析中使用到的小波函數(shù)具有不唯一性，即小波函數(shù) 具有多樣性。小波分析在工程應(yīng)用中，一個(gè)十分重要的問題就是最優(yōu)小波基的選擇問題，因?yàn)橛貌煌男〔ɑ治鐾粋€(gè)問題會(huì)產(chǎn)生不同的結(jié)果。目前我們主要是通過用小波分析方法處理信號(hào)的結(jié)果與理論結(jié)果的誤差來判定小波基的好壞，由此決定小波基。
???? 常用小波基有Haar小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet小波、Meyer小波等。
?????1.Haar小波Haar函數(shù)是小波分析中最早用到的一個(gè)具有緊支撐的正交小波函數(shù)，也是最簡單的一個(gè)小波函數(shù)，它是支撐域在范圍內(nèi)的單個(gè)矩形波。Haar函數(shù)的定義如下：Haar小波在時(shí)域上是不連續(xù)的，所以作為基本小波性能不是特別好。但它也有自己的優(yōu)點(diǎn):計(jì)算簡單。不但與正交，而且與自己的整數(shù)位移正交，因此，在的多分辨率系統(tǒng)中，Haar小波構(gòu)成一組最簡單的正交歸一的小波族。
?[phi,g1,xval] = wavefun('haar',20);
subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t')；
title('haar 時(shí)域');
g2=fft(g1);
g3=abs(g2);
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f')title('haar 頻域')；
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
2.Daubechies(dbN)小波Daubechies小波是世界著名的小波分析學(xué)者Inrid·Daubechies構(gòu)造的小波函數(shù)，簡寫為dbN，N是小波的階數(shù)。小波和尺度函數(shù)中的支撐區(qū)為的消失矩為。除（Harr小波）外，dbN不具有對(duì)稱性（即非線性相位）。除（Harr小波）外，dbN沒有明確的表達(dá)式，但轉(zhuǎn)換函數(shù)h的平方模是明確的:令，其中為二項(xiàng)式的系數(shù)，則有其中：Daubechies小波具有以下特點(diǎn)：在時(shí)域是有限支撐的，即長度有限。在頻域在處有N階零點(diǎn)。和它的整數(shù)位移正交歸一，即。小波函數(shù)可以由所謂“尺度函數(shù)”求出來。尺度函數(shù)為低通函數(shù)，長度有限，支撐域在的范圍內(nèi)。
?db4的時(shí)域和頻域波形：
[phi,g1,xval] = wavefun('db4',10);
subplot(2,1,1);
plot(xval,g1,'LineWidth',2);
xlabel('t')title('db4 時(shí)域');
g2=fft(g1);
g3=abs(g2);
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f')title('db4 頻域')；
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
?Daubechies小波常用來分解和重構(gòu)信號(hào)，作為濾波器使用:
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db4'); %計(jì)算該小波的4個(gè)濾波器
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('重構(gòu)低通濾波器');
subplot(2,2,4);
stem(Hi_R,'LineWidth',2);
title('重構(gòu)高通濾波器');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
? 3.Mexican Hat(mexh)小波Mexican Hat函數(shù)為Gauss函數(shù)的二階導(dǎo)數(shù)：因?yàn)樗男螤钕衲鞲缑钡慕孛?#xff0c;所以也稱為墨西哥帽函數(shù)。Mexihat小波的時(shí)域和頻域波形：
?Mexihat小波的時(shí)域和頻域波形：
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 時(shí)域');
g2=fft(g1);
g3=(abs(g2));
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 頻域');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
Mexihat小波的特點(diǎn)：在時(shí)間域與頻率域都有很好的局部化，并且滿足。不存在尺度函數(shù)，所以Mexihat小波函數(shù)不具有正交性。
?4.Morlet小波它是高斯包絡(luò)下的單頻率副正弦函數(shù)：其中C是重構(gòu)時(shí)的歸一化常數(shù)。Morlet小波沒有尺度函數(shù)，而且是非正交分解。Morlet小波的時(shí)域和頻域波形圖：
?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 時(shí)域');
g2=fft(g1);
g3=(abs(g2));
subplot(2,1,2);
plot(g3,'LineWidth',2);
xlabel('f');
title('mexihat 頻域');
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
注:
獲取小波系數(shù)的兩個(gè)函數(shù)理解：
Wpcoef?是求解某個(gè)節(jié)點(diǎn)的小波包系數(shù)，數(shù)據(jù)長度是1/(2^n)（分解n層的話），其實(shí)就是將原始信號(hào)化成2^N段，每段的長度是相等的且比原信號(hào)短
wprcoef是把某個(gè)節(jié)點(diǎn)的小波包系數(shù)重構(gòu)，得到的是和原信號(hào)一樣長度的信號(hào)。
總結(jié)
以上是生活随笔為你收集整理的小波与小波包、小波包分解与信号重构、小波包能量特征提取 暨 小波包分解后实现按频率大小分布重新排列(Matlab 程序详解）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： PS学习笔记--去掉图片上不想要的部分
- 下一篇： 最新去广告扒站仿站网源码