當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

小波包分解与重构

發(fā)布時(shí)間：2023/12/31 编程问答 28 豆豆

生活随笔收集整理的這篇文章主要介紹了小波包分解与重构小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

1、小波變換的理解

傅里葉變換——短時(shí)傅里葉變換——小波變換。

參考文獻(xiàn)：以下兩篇參考資料講述得十分清楚，有助于理解小波變換。

但具體的數(shù)學(xué)角度闡述，請(qǐng)參考其他資料。

（1）知乎專欄：形象易懂講解算法I——小波變換

https://zhuanlan.zhihu.com/p/22450818

（2）知乎專欄：傅里葉分析之掐死教程。

https://zhuanlan.zhihu.com/p/19763358

2、小波包分解

小波包是為了克服小波分解在高頻段的頻率分辨率較差，而在低頻段的時(shí)間分辨率較差的問題的基礎(chǔ)上而提出的。

它是一種更精細(xì)的信號(hào)分析的方法，提高了信號(hào)的時(shí)域分辨率。

下面是兩者的對(duì)比圖：

3、能量譜

??????基于小波包分解提取多尺度空間能量特征的原理是把不同分解尺度上的信號(hào)能量求解出來，將這些能量值按尺度順序排列成特征向量供識(shí)別使用。

20180510補(bǔ)充更新：具體計(jì)算公式如下所示，本文中未使用重構(gòu)后的系數(shù)進(jìn)行能量值計(jì)算，直接使用小波包分解后的系數(shù)，參考文獻(xiàn)《基于小波包能量特征的滾動(dòng)軸承故障監(jiān)測(cè)方法》。

4、Matlab代碼

給出兩部分代碼，寫成兩個(gè)函數(shù)。一個(gè)是小波包分解與重構(gòu)，另一個(gè)是能量譜函數(shù)。

下載地址：https://download.csdn.net/download/ckzhb/10030651

代碼名稱：wavelet_packetdecomposition_reconstruct

function wpt= wavelet_packetdecomposition_reconstruct( x,n,wpname )
%% 對(duì)信號(hào)進(jìn)行小波包分解，得到節(jié)點(diǎn)的小波包系數(shù)。然后對(duì)每個(gè)節(jié)點(diǎn)系數(shù)進(jìn)行重構(gòu)。?
% 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('原始信號(hào)');
subplot(5,1,2);
plot(cfs0);
title(['結(jié)點(diǎn) ',num2str(n) ' ?1',' 系數(shù)'])
subplot(5,1,3);
plot(cfs1);
title(['結(jié)點(diǎn) ',num2str(n) ' ?2',' 系數(shù)'])
subplot(5,1,4);
plot(cfs2);
title(['結(jié)點(diǎn) ',num2str(n) ' ?3',' 系數(shù)'])
subplot(5,1,5);
plot(cfs3);
title(['結(jié)點(diǎn) ',num2str(n) ' ?4',' 系數(shù)'])
%% 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('原始信號(hào)');
subplot(5,1,2);
plot(rex0);
title(['重構(gòu)結(jié)點(diǎn) ',num2str(n) ' ?1',' 系數(shù)'])
subplot(5,1,3);
plot(rex1);
title(['重構(gòu)結(jié)點(diǎn) ',num2str(n) ' ?2',' 系數(shù)'])
subplot(5,1,4);
plot(rex2);
title(['重構(gòu)結(jié)點(diǎn) ',num2str(n) ' ?3',' 系數(shù)'])
subplot(5,1,5);
plot(rex3);
title(['重構(gòu)結(jié)點(diǎn) ',num2str(n) ' ?4',' 系數(shù)'])
end

代碼名稱：wavelet_energy_spectrum
function E = wavelet_energy_spectrum( wpt,n )
%% 計(jì)算每一層每一個(gè)節(jié)點(diǎn)的能量
% ?wpt-wavelet packet tree
% ?n-第n層能量
%?
% Author hubery_zhang
% Date ?20170714

%%
% 求第n層第i個(gè)節(jié)點(diǎn)的系數(shù)
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
%求每個(gè)節(jié)點(diǎn)的概率
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('結(jié)點(diǎn)');
ylabel('能量百分比/%');
for j=1:2^n
text(x(j),p_node(i),num2str(p_node(j),'%0.2f'),...
? ? 'HorizontalAlignment','center',...
? ? 'VerticalAlignment','bottom')
end

end
————————————————
版權(quán)聲明：本文為CSDN博主「hubery_zhang」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議，轉(zhuǎn)載請(qǐng)附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/ckzhb/article/details/78288847

?小波與小波包、小波包分解與信號(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ù)