热门标签 | HotTags
当前位置:  开发笔记 > 编程语言 > 正文

【转】小波与小波包、小波包分解与信号重构、小波包能量特征提取暨小波包分解后实现按频率大小分布重新排列(Matlab程序详解)

转:https:blog.csdn.netcqfdcwarticledetails84995904    小波与小波包、小波包分解与信号重构、小波包能量特征提取  (Matlab程

转:https://blog.csdn.net/cqfdcw/article/details/84995904


        小波与小波包、小波包分解与信号重构、小波包能量特征提取   (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的正弦波混合的,我们用小波包来分解。




  1.  
    clear all


  2.  
    clc


  3.  
    fs=1024; %采样频率


  4.  
    f1=100; %信号的第一个频率


  5.  
    f2=300; %信号第二个频率


  6.  
    t=0:1/fs:1;


  7.  
    s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %生成混合信号


  8.  
    [tt]=wpdec(s,3,'dmey'); %小波包分解,3代表分解3层,'dmey'使用meyr小波


  9.  
    plot(tt) %画小波包树图


  10.  
    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进来就行了。




  1.  
    t=wpdec(wave,3,'dmey');


  2.  
    t2 = wpjoin(t,[3;4;5;6]);


  3.  
    sNod = read(t,'sizes',[3,4,5,6]);


  4.  
    cfs3 = zeros(sNod(1,:));


  5.  
    cfs4 = zeros(sNod(2,:));


  6.  
    cfs5 = zeros(sNod(3,:));


  7.  
    cfs6 = zeros(sNod(4,:));


  8.  
    t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);


  9.  
    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




  1.  
    t=wpdec(wave,3,'dmey');


  2.  
    t2 = wpjoin(t,[3;4;5;6]);


  3.  
    cfs3=wpcoef(t,3);


  4.  
    cfs4=wpcoef(t,4);


  5.  
    cfs5=wpcoef(t,5);


  6.  
    cfs6=wpcoef(t,6);


  7.  
    t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);


  8.  
    wave2=wprec(t3);

解释:

第一行:将wave 用 dmey小波进行3层小波包分解,获得一个小波包树 t

第二行:将小波包树的第二行的四个节点收起来,也就是让第二行的节点变为树的最底层节点。

第三~六行:获取四个节点的小波包系数 (小波包系数就是一个一维向量)

第七行:将四个节点的系数重组到t3小波树中

第八行:对t3小波树进行重构,获得信号wave2

可以看出,该例子就是对一个小波包展开了,又原封不动的装回去了,所以说 wave2和wave是一样的。

注意,wpjoin命令在这里是必要的,因为write函数只能将最底层的节点写进去。也就是说,如果我们将第三层的小波包系数进行修改的话,就不用wpjoin了,具体可以看例3

例3




  1.  
    t=wpdec(wave,3,'dmey');


  2.  
    cfs7=wpcoef(t,7);


  3.  
    cfs8=wpcoef(t,8);


  4.  
    cfs9=wpcoef(t,9);


  5.  
    cfs10=wpcoef(t,10);


  6.  
    cfs11=wpcoef(t,11);


  7.  
    cfs12=wpcoef(t,12);


  8.  
    cfs13=wpcoef(t,13);


  9.  
    cfs14=wpcoef(t,14);


  10.  
    t3=write(t,'cfs',7,cfs7,'cfs',8,cfs8,'cfs',9,cfs9,'cfs',10,cfs10,'cfs',11,cfs11,'cfs',...


  11.  
    12,cfs12,'cfs',13,cfs13,'cfs',14,cfs14);


  12.  
    y=wprec(t3);

该例子也是对一个小波包展开了,又原封不动的装回去了,只不过这次是直接对第三层节点进行的。


5.小波包-----信号分解与重构(方法2)

该方法只能对某一节点信号系数分别进行重构,不能实现多个节点系数组合进行重构.

接下来进行举例说明:




  1.  
    x_input=x_train(:,1,1); %输入数据


  2.  
    plot(x_input);title('输入信号时域图像') %绘制输入信号时域图像




  1.  
    %% 查看频谱范围


  2.  
    x=x_input;


  3.  
    fs=128;


  4.  
    N=length(x); %采样点个数


  5.  
    signalFFT=abs(fft(x,N));%真实的幅值


  6.  
    Y=2*signalFFT/N;


  7.  
    f=(0:N/2)*(fs/N);


  8.  
    figure;plot(f,Y(1:N/2+1));


  9.  
    ylabel('amp'); xlabel('frequency');title('输入信号的频谱');grid on

接下来进行4层小波包分解




  1.  
    wpt=wpdec(x_input,3,'dmey'); %进行3层小波包分解


  2.  
    plot(wpt); %绘制小波包树

接下来,我们查看第3层8个节点的频谱分布(我的输入信号采样频率是128,按照采样定理小波包分解根节点(0,0)处的频率应该为0-64Hz,按照这个计算(3,0)节点为0-8Hz,后面依次以8Hz为一个段递增)

首先展示最初的重构方法(频率排布顺序混乱,常见理解错误)




  1.  
    for i=0:7


  2.  
    rex3(:,i+1)=wprcoef(wpt,[3 i]); %实现对节点小波节点进行重构


  3.  
    end


  4.  
     


  5.  
    figure; %绘制第3层各个节点分别重构后信号的频谱


  6.  
    for i=0:7


  7.  
    subplot(2,4,i+1);


  8.  
    x_sign=rex3(:,i+1);


  9.  
    N=length(x_sign); %采样点个数


  10.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值


  11.  
    Y=2*signalFFT/N;


  12.  
    f=(0:N/2)*(fs/N);


  13.  
    plot(f,Y(1:N/2+1));


  14.  
    ylabel('amp'); xlabel('frequency');grid on


  15.  
    axis([0 50 0 0.03]); title(['小波包第3层',num2str(i),'节点信号频谱']);


  16.  
    end

绘制完后你会发现频谱分布并不是按照之前理解的顺序依次递增排列,如下所示

那么问题出在哪里呢?经过仔细深入验证后发现问题出在小波包节点的频谱划分“不是严格按照上述理解的顺序排列的”(可能是一种格雷编码或者其他),要解决这个问题我们需要对节点顺序进行重新编码排序。参考这里https://www.ilovematlab.cn/thread-122226-1-1.html?tdsourcetag=s_pctim_aiomsg




  1.  
    nodes=[7;8;9;10;11;12;13;14]; %第3层的节点号


  2.  
    ord=wpfrqord(nodes);  %小波包系数重排,ord是重排后小波包系数索引构成的矩阵 如3层分解的[1;2;4;3;7;8;6;5]


  3.  
    nodes_ord=nodes(ord); %重排后的小波系数

注意:节点的命名规则是从(1,0)开始,叫1号, (1,1)是2号………依此类推,(3,0)是7号,(3,7)是14号。

那么我们来看看经过上面这段重排序运行后nodes_ord中顺序是什么?

接下来我们再绘制一下第三层8个节点重构信号的频谱如下




  1.  
    for i=1:8


  2.  
    rex3(:,i)=wprcoef(wpt,nodes_ord(i)); %实现对节点小波节点进行重构


  3.  
    end


  4.  
     


  5.  
    figure; %绘制第3层各个节点分别重构后信号的频谱


  6.  
    for i=0:7


  7.  
    subplot(2,4,i+1);


  8.  
    x_sign= rex3(:,i+1);


  9.  
    N=length(x_sign); %采样点个数


  10.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值


  11.  
    Y=2*signalFFT/N;


  12.  
    f=(0:N/2)*(fs/N);


  13.  
    plot(f,Y(1:N/2+1));


  14.  
    ylabel('amp'); xlabel('frequency');grid on


  15.  
    axis([0 50 0 0.03]); title(['小波包第3层',num2str(i),'节点信号频谱']);


  16.  
    end

到这里为止不知道大家有没有明白我要表达的意思,如果没有明白可以再反复看理解,或者自己进行分解后观察每个节点的频谱后或许就理解了。


6.小波包分解------能量特征提取(方法1)

接下来绘制第3层各个频段的能量占比




  1.  
    %% wavelet packet coefficients. 求取小波包分解的各个节点的小波包系数


  2.  
    cfs3_0=wpcoef(wpt,nodes_ord(1)); %对重排序后第3层0节点的小波包系数0-8Hz


  3.  
    cfs3_1=wpcoef(wpt,nodes_ord(2)); %对重排序后第3层0节点的小波包系数8-16Hz


  4.  
    cfs3_2=wpcoef(wpt,nodes_ord(3)); %对重排序后第3层0节点的小波包系数16-24Hz


  5.  
    cfs3_3=wpcoef(wpt,nodes_ord(4)); %对重排序后第3层0节点的小波包系数24-32Hz


  6.  
    cfs3_4=wpcoef(wpt,nodes_ord(5)); %对重排序后第3层0节点的小波包系数32-40Hz


  7.  
    cfs3_5=wpcoef(wpt,nodes_ord(6)); %对重排序后第3层0节点的小波包系数40-48Hz


  8.  
    cfs3_6=wpcoef(wpt,nodes_ord(7)); %对重排序后第3层0节点的小波包系数48-56Hz


  9.  
    cfs3_7=wpcoef(wpt,nodes_ord(8)); %对重排序后第3层0节点的小波包系数56-64Hz


  10.  
     


  11.  
    E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范数:就是norm(...,1),即各元素绝对值之和;2-范数:就是norm(...,2),即各元素平方和开根号;


  12.  
    E_cfs3_1=norm(cfs3_1,2)^2;


  13.  
    E_cfs3_2=norm(cfs3_2,2)^2;


  14.  
    E_cfs3_3=norm(cfs3_3,2)^2;


  15.  
    E_cfs3_4=norm(cfs3_4,2)^2;


  16.  
    E_cfs3_5=norm(cfs3_5,2)^2;


  17.  
    E_cfs3_6=norm(cfs3_6,2)^2;


  18.  
    E_cfs3_7=norm(cfs3_7,2)^2;


  19.  
    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;


  20.  
     


  21.  
    p_node(0)= 100*E_cfs3_0/E_total; % 求得每个节点的占比


  22.  
    p_node(2)= 100*E_cfs3_1/E_total; % 求得每个节点的占比


  23.  
    p_node(3)= 100*E_cfs3_2/E_total; % 求得每个节点的占比


  24.  
    p_node(4)= 100*E_cfs3_3/E_total; % 求得每个节点的占比


  25.  
    p_node(5)= 100*E_cfs3_4/E_total; % 求得每个节点的占比


  26.  
    p_node(6)= 100*E_cfs3_5/E_total; % 求得每个节点的占比


  27.  
    p_node(7)= 100*E_cfs3_6/E_total; % 求得每个节点的占比


  28.  
    p_node(8)= 100*E_cfs3_7/E_total; % 求得每个节点的占比


  29.  
     


  30.  
    figure;


  31.  
    x=1:8;


  32.  
    bar(x,p_node);


  33.  
    title('各个频段能量所占的比例');


  34.  
    xlabel('频率 Hz');


  35.  
    ylabel('能量百分比/%');


  36.  
    for j=1:8


  37.  
    text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...


  38.  
    'HorizontalAlignment','center',...


  39.  
    'VerticalAlignment','bottom')


  40.  
    end


7.小波包分解------能量特征提取(方法2)

直接运行matlab自带函数,如下

E = wenergy(wpt);   %该函数只能对最后一层(底层)节点进行能量提取

(引用:https://blog.csdn.net/it_beecoder/article/details/78668273)


这里对比一下,和方法1得到的图形基本上是一致的,不同之处在于排列顺序变了,方法2中的顺序是按照小波包分解函数自动排列顺序,方法1经过重排序后为按照频率段递增顺序依次排列顺序的




  1.  
    %% 绘制重构前各个频段小波包系数


  2.  
    figure(1);


  3.  
    subplot(4,1,1);


  4.  
    plot(x_input);


  5.  
    title('原始信号');


  6.  
    subplot(4,1,2);


  7.  
    plot(cfs3_0);


  8.  
    title(['层数 ',num2str(3) ' 节点 0的小波0-8Hz',' 系数'])


  9.  
    subplot(4,1,3);


  10.  
    plot(cfs3_1);


  11.  
    title(['层数 ',num2str(3) ' 节点 1的小波8-16Hz',' 系数'])


  12.  
    subplot(4,1,4);


  13.  
    plot(cfs3_2);


  14.  
    title(['层数 ',num2str(3) ' 节点 2的小波16-24Hz',' 系数'])




  1.  
    %% 绘制重构后时域各个特征频段的图形


  2.  
    figure(3);


  3.  
    subplot(3,1,1);


  4.  
    plot(rex3(:,1));


  5.  
    title('重构后0-8Hz频段信号');


  6.  
    subplot(3,1,2);


  7.  
    plot(rex3(:,2));


  8.  
    title('重构后8-16Hz频段信号')


  9.  
    subplot(3,1,3);


  10.  
    plot(rex3(:,3));


  11.  
    title('重构后16-24Hz频段信号');

以上过程完整代码




  1.  
    clear all;


  2.  
     


  3.  
    x_input=x_train(:,1,1); %输入数据


  4.  
    plot(x_input);title('输入信号时域图像') %绘制输入信号时域图像


  5.  
     


  6.  
    x=x_input; % 查看频谱范围


  7.  
    fs=128;


  8.  
    N=length(x); %采样点个数


  9.  
    signalFFT=abs(fft(x,N));%真实的幅值


  10.  
    Y=2*signalFFT/N;


  11.  
    f=(0:N/2)*(fs/N);


  12.  
    figure;plot(f,Y(1:N/2+1));


  13.  
    ylabel('amp'); xlabel('frequency');title('输入信号的频谱');grid on


  14.  
     


  15.  
    wpt=wpdec(x_input,3,'dmey'); %进行4层小波包分解


  16.  
    plot(wpt); %绘制小波包树


  17.  
     


  18.  
    %% 实现对节点顺序按照频率递增进行重排序


  19.  
    % nodes=get(wpt,'tn'); %小波包分解系数 为什么nodes是[7;8;9;10;11;12;13;14]


  20.  
    % N_cfs=length(nodes); %小波包系数个数


  21.  
    nodes=[7;8;9;10;11;12;13;14];


  22.  
    ord=wpfrqord(nodes); %小波包系数重排,ord是重排后小波包系数索引构成的矩阵 如3层分解的[1;2;4;3;7;8;6;5]


  23.  
    nodes_ord=nodes(ord); %重排后的小波系数


  24.  
     


  25.  
    %% 实现对节点小波节点进行重构


  26.  
    for i=1:8


  27.  
    rex3(:,i)=wprcoef(wpt,nodes_ord(i));


  28.  
    end


  29.  
     


  30.  
    %% 绘制第3层各个节点分别重构后信号的频谱


  31.  
    figure;


  32.  
    for i=0:7


  33.  
    subplot(2,4,i+1);


  34.  
    x_sign= rex3(:,i+1);


  35.  
    N=length(x_sign); %采样点个数


  36.  
    signalFFT=abs(fft(x_sign,N));%真实的幅值


  37.  
    Y=2*signalFFT/N;


  38.  
    f=(0:N/2)*(fs/N);


  39.  
    plot(f,Y(1:N/2+1));


  40.  
    ylabel('amp'); xlabel('frequency');grid on


  41.  
    axis([0 50 0 0.03]); title(['小波包第3层',num2str(i),'节点信号频谱']);


  42.  
    end


  43.  
     


  44.  
    %% wavelet packet coefficients. 求取小波包分解的各个频段的小波包系数


  45.  
    cfs3_0=wpcoef(wpt,nodes_ord(1)); %对重排序后第3层0节点的小波包系数0-8Hz


  46.  
    cfs3_1=wpcoef(wpt,nodes_ord(2)); %对重排序后第3层0节点的小波包系数8-16Hz


  47.  
    cfs3_2=wpcoef(wpt,nodes_ord(3)); %对重排序后第3层0节点的小波包系数16-24Hz


  48.  
    cfs3_3=wpcoef(wpt,nodes_ord(4)); %对重排序后第3层0节点的小波包系数24-32Hz


  49.  
    cfs3_4=wpcoef(wpt,nodes_ord(5)); %对重排序后第3层0节点的小波包系数32-40Hz


  50.  
    cfs3_5=wpcoef(wpt,nodes_ord(6)); %对重排序后第3层0节点的小波包系数40-48Hz


  51.  
    cfs3_6=wpcoef(wpt,nodes_ord(7)); %对重排序后第3层0节点的小波包系数48-56Hz


  52.  
    cfs3_7=wpcoef(wpt,nodes_ord(8)); %对重排序后第3层0节点的小波包系数56-64Hz


  53.  
    %% 求取小波包分解的各个频段的能量


  54.  
    E_cfs3_0=norm(cfs3_0,2)^2; %% 1-范数:就是norm(...,1),即各元素绝对值之和;2-范数:就是norm(...,2),即各元素平方和开根号;


  55.  
    E_cfs3_1=norm(cfs3_1,2)^2;


  56.  
    E_cfs3_2=norm(cfs3_2,2)^2;


  57.  
    E_cfs3_3=norm(cfs3_3,2)^2;


  58.  
    E_cfs3_4=norm(cfs3_4,2)^2;


  59.  
    E_cfs3_5=norm(cfs3_5,2)^2;


  60.  
    E_cfs3_6=norm(cfs3_6,2)^2;


  61.  
    E_cfs3_7=norm(cfs3_7,2)^2;


  62.  
    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;


  63.  
     


  64.  
    p_node(0)= 100*E_cfs3_0/E_total; % 求得每个节点的占比


  65.  
    p_node(2)= 100*E_cfs3_1/E_total; % 求得每个节点的占比


  66.  
    p_node(3)= 100*E_cfs3_2/E_total; % 求得每个节点的占比


  67.  
    p_node(4)= 100*E_cfs3_3/E_total; % 求得每个节点的占比


  68.  
    p_node(5)= 100*E_cfs3_4/E_total; % 求得每个节点的占比


  69.  
    p_node(6)= 100*E_cfs3_5/E_total; % 求得每个节点的占比


  70.  
    p_node(7)= 100*E_cfs3_6/E_total; % 求得每个节点的占比


  71.  
    p_node(8)= 100*E_cfs3_7/E_total; % 求得每个节点的占比


  72.  
     


  73.  
    figure;


  74.  
    x=1:8;


  75.  
    bar(x,p_node);


  76.  
    title('各个频段能量所占的比例');


  77.  
    xlabel('频率 Hz');


  78.  
    ylabel('能量百分比/%');


  79.  
    for j=1:8


  80.  
    text(x(j),p_node(j),num2str(p_node(j),'%0.2f'),...


  81.  
    'HorizontalAlignment','center',...


  82.  
    'VerticalAlignment','bottom')


  83.  
    end


  84.  
    % E = wenergy(wpt); %求取各个节点能量


  85.  
     


  86.  
    %% 绘制重构前各个特征频段小波包系数的图形


  87.  
    figure(1);


  88.  
    subplot(4,1,1);


  89.  
    plot(x_input);


  90.  
    title('原始信号');


  91.  
    subplot(4,1,2);


  92.  
    plot(cfs3_0);


  93.  
    title(['层数 ',num2str(3) ' 节点 0的小波0-8Hz',' 系数'])


  94.  
    subplot(4,1,3);


  95.  
    plot(cfs3_1);


  96.  
    title(['层数 ',num2str(3) ' 节点 1的小波8-16Hz',' 系数'])


  97.  
    subplot(4,1,4);


  98.  
    plot(cfs3_2);


  99.  
    title(['层数 ',num2str(3) ' 节点 2的小波16-24Hz',' 系数'])


  100.  
     


  101.  
    %% 绘制重构后时域各个特征频段的图形


  102.  
    figure(3);


  103.  
    subplot(3,1,1);


  104.  
    plot(rex3(:,1));


  105.  
    title('重构后0-8Hz频段信号');


  106.  
    subplot(3,1,2);


  107.  
    plot(rex3(:,2));


  108.  
    title('重构后8-16Hz频段信号')


  109.  
    subplot(3,1,3);


  110.  
    plot(rex3(:,3));


  111.  
    title('重构后16-24Hz频段信号');


  112.  
     


8.小波----常见基函数

     与标准的傅里叶变换相比,小波分析中使用到的小波函数具有不唯一性,即小波函数 具有多样性。小波分析在工程应用中,一个十分重要的问题就是最优小波基的选择问题,因为用不同的小波基分析同一个问题会产生不同的结果。目前我们主要是通过用小波分析方法处理信号的结果与理论结果的误差来判定小波基的好坏,由此决定小波基。

     常用小波基有Haar小波、Daubechies(dbN)小波、Mexican Hat(mexh)小波、Morlet小波、Meyer小波等。

     1.Haar小波Haar函数是小波分析中最早用到的一个具有紧支撑的正交小波函数,也是最简单的一个小波函数,它是支撑域在范围内的单个矩形波。Haar函数的定义如下:Haar小波在时域上是不连续的,所以作为基本小波性能不是特别好。但它也有自己的优点:计算简单。不但与正交,而且与自己的整数位移正交,因此,在的多分辨率系统中,Haar小波构成一组最简单的正交归一的小波族。




  1.  
    [phi,g1,xval] = wavefun('haar',20);


  2.  
    subplot(2,1,1);


  3.  
    plot(xval,g1,'LineWidth',2);


  4.  
    xlabel('t');


  5.  
    title('haar 时域');


  6.  
    g2=fft(g1);


  7.  
    g3=abs(g2);


  8.  
    subplot(2,1,2);


  9.  
    plot(g3,'LineWidth',2);


  10.  
    xlabel('f')title('haar 频域');

                                      

2.Daubechies(dbN)小波Daubechies小波是世界著名的小波分析学者Inrid·Daubechies构造的小波函数,简写为dbN,N是小波的阶数。小波和尺度函数中的支撑区为的消失矩为。除(Harr小波)外,dbN不具有对称性(即非线性相位)。除(Harr小波)外,dbN没有明确的表达式,但转换函数h的平方模是明确的:令,其中为二项式的系数,则有其中:Daubechies小波具有以下特点:在时域是有限支撑的,即长度有限。在频域在处有N阶零点。和它的整数位移正交归一,即。小波函数可以由所谓“尺度函数”求出来。尺度函数为低通函数,长度有限,支撑域在的范围内。




  1.  
    db4的时域和频域波形:


  2.  
    [phi,g1,xval] = wavefun('db4',10);


  3.  
    subplot(2,1,1);


  4.  
    plot(xval,g1,'LineWidth',2);


  5.  
    xlabel('t')title('db4 时域');


  6.  
    g2=fft(g1);


  7.  
    g3=abs(g2);


  8.  
    subplot(2,1,2);


  9.  
    plot(g3,'LineWidth',2);


  10.  
    xlabel('f')title('db4 频域');

                                  




  1.  
    Daubechies小波常用来分解和重构信号,作为滤波器使用:


  2.  
    [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db4'); %计算该小波的4个滤波器


  3.  
    subplot(2,2,1);


  4.  
    stem(Lo_D,'LineWidth',2);


  5.  
    title('分解低通滤波器');


  6.  
    subplot(2,2,2);


  7.  
    stem(Hi_D,'LineWidth',2);


  8.  
    title('分解高通滤波器');


  9.  
    subplot(2,2,3);


  10.  
    stem(Lo_R,'LineWidth',2);


  11.  
    title('重构低通滤波器');


  12.  
    subplot(2,2,4);


  13.  
    stem(Hi_R,'LineWidth',2);


  14.  
    title('重构高通滤波器');

                              

  3.Mexican Hat(mexh)小波Mexican Hat函数为Gauss函数的二阶导数:因为它的形状像墨西哥帽的截面,所以也称为墨西哥帽函数。Mexihat小波的时域和频域波形:




  1.  
    Mexihat小波的时域和频域波形:


  2.  
    d=-6; h=6; n=100;


  3.  
    [g1,x]=mexihat(d,h,n);


  4.  
    subplot(2,1,1);


  5.  
    plot(x,g1,'LineWidth',2);


  6.  
    xlabel('t');


  7.  
    title('Mexihat 时域');


  8.  
    g2=fft(g1);


  9.  
    g3=(abs(g2));


  10.  
    subplot(2,1,2);


  11.  
    plot(g3,'LineWidth',2);


  12.  
    xlabel('f');


  13.  
    title('mexihat 频域');

                                  

Mexihat小波的特点:在时间域与频率域都有很好的局部化,并且满足。不存在尺度函数,所以Mexihat小波函数不具有正交性。

 4.Morlet小波它是高斯包络下的单频率副正弦函数:其中C是重构时的归一化常数。Morlet小波没有尺度函数,而且是非正交分解。Morlet小波的时域和频域波形图:




  1.  
    d=-6; h=6; n=100;


  2.  
    [g1,x]=morlet(d,h,n);


  3.  
    subplot(2,1,1);


  4.  
    plot(x,g1,'LineWidth',2);


  5.  
    xlabel('t');


  6.  
    title('morlet 时域');


  7.  
    g2=fft(g1);


  8.  
    g3=(abs(g2));


  9.  
    subplot(2,1,2);


  10.  
    plot(g3,'LineWidth',2);


  11.  
    xlabel('f');


  12.  
    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



推荐阅读
  • 本文介绍了九度OnlineJudge中的1002题目“Grading”的解决方法。该题目要求设计一个公平的评分过程,将每个考题分配给3个独立的专家,如果他们的评分不一致,则需要请一位裁判做出最终决定。文章详细描述了评分规则,并给出了解决该问题的程序。 ... [详细]
  • 本文介绍了P1651题目的描述和要求,以及计算能搭建的塔的最大高度的方法。通过动态规划和状压技术,将问题转化为求解差值的问题,并定义了相应的状态。最终得出了计算最大高度的解法。 ... [详细]
  • 微软头条实习生分享深度学习自学指南
    本文介绍了一位微软头条实习生自学深度学习的经验分享,包括学习资源推荐、重要基础知识的学习要点等。作者强调了学好Python和数学基础的重要性,并提供了一些建议。 ... [详细]
  • 这是原文链接:sendingformdata许多情况下,我们使用表单发送数据到服务器。服务器处理数据并返回响应给用户。这看起来很简单,但是 ... [详细]
  • 本文主要解析了Open judge C16H问题中涉及到的Magical Balls的快速幂和逆元算法,并给出了问题的解析和解决方法。详细介绍了问题的背景和规则,并给出了相应的算法解析和实现步骤。通过本文的解析,读者可以更好地理解和解决Open judge C16H问题中的Magical Balls部分。 ... [详细]
  • 推荐系统遇上深度学习(十七)详解推荐系统中的常用评测指标
    原创:石晓文小小挖掘机2018-06-18笔者是一个痴迷于挖掘数据中的价值的学习人,希望在平日的工作学习中,挖掘数据的价值, ... [详细]
  • 解决VS写C#项目导入MySQL数据源报错“You have a usable connection already”问题的正确方法
    本文介绍了在VS写C#项目导入MySQL数据源时出现报错“You have a usable connection already”的问题,并给出了正确的解决方法。详细描述了问题的出现情况和报错信息,并提供了解决该问题的步骤和注意事项。 ... [详细]
  • 本文详细介绍了Java中vector的使用方法和相关知识,包括vector类的功能、构造方法和使用注意事项。通过使用vector类,可以方便地实现动态数组的功能,并且可以随意插入不同类型的对象,进行查找、插入和删除操作。这篇文章对于需要频繁进行查找、插入和删除操作的情况下,使用vector类是一个很好的选择。 ... [详细]
  • 本文详细介绍了MySQL表分区的创建、增加和删除方法,包括查看分区数据量和全库数据量的方法。欢迎大家阅读并给予点评。 ... [详细]
  • 猜字母游戏
    猜字母游戏猜字母游戏——设计数据结构猜字母游戏——设计程序结构猜字母游戏——实现字母生成方法猜字母游戏——实现字母检测方法猜字母游戏——实现主方法1猜字母游戏——设计数据结构1.1 ... [详细]
  • [大整数乘法] java代码实现
    本文介绍了使用java代码实现大整数乘法的过程,同时也涉及到大整数加法和大整数减法的计算方法。通过分治算法来提高计算效率,并对算法的时间复杂度进行了研究。详细代码实现请参考文章链接。 ... [详细]
  • 本文介绍了一些Java开发项目管理工具及其配置教程,包括团队协同工具worktil,版本管理工具GitLab,自动化构建工具Jenkins,项目管理工具Maven和Maven私服Nexus,以及Mybatis的安装和代码自动生成工具。提供了相关链接供读者参考。 ... [详细]
  • 【shell】网络处理:判断IP是否在网段、两个ip是否同网段、IP地址范围、网段包含关系
    本文介绍了使用shell脚本判断IP是否在同一网段、判断IP地址是否在某个范围内、计算IP地址范围、判断网段之间的包含关系的方法和原理。通过对IP和掩码进行与计算,可以判断两个IP是否在同一网段。同时,还提供了一段用于验证IP地址的正则表达式和判断特殊IP地址的方法。 ... [详细]
  • 模板引擎StringTemplate的使用方法和特点
    本文介绍了模板引擎StringTemplate的使用方法和特点,包括强制Model和View的分离、Lazy-Evaluation、Recursive enable等。同时,还介绍了StringTemplate语法中的属性和普通字符的使用方法,并提供了向模板填充属性的示例代码。 ... [详细]
  • 阿里Treebased Deep Match(TDM) 学习笔记及技术发展回顾
    本文介绍了阿里Treebased Deep Match(TDM)的学习笔记,同时回顾了工业界技术发展的几代演进。从基于统计的启发式规则方法到基于内积模型的向量检索方法,再到引入复杂深度学习模型的下一代匹配技术。文章详细解释了基于统计的启发式规则方法和基于内积模型的向量检索方法的原理和应用,并介绍了TDM的背景和优势。最后,文章提到了向量距离和基于向量聚类的索引结构对于加速匹配效率的作用。本文对于理解TDM的学习过程和了解匹配技术的发展具有重要意义。 ... [详细]
author-avatar
泽旺多吉外_680
这个家伙很懒,什么也没留下!
PHP1.CN | 中国最专业的PHP中文社区 | DevBox开发工具箱 | json解析格式化 |PHP资讯 | PHP教程 | 数据库技术 | 服务器技术 | 前端开发技术 | PHP框架 | 开发工具 | 在线工具
Copyright © 1998 - 2020 PHP1.CN. All Rights Reserved | 京公网安备 11010802041100号 | 京ICP备19059560号-4 | PHP1.CN 第一PHP社区 版权所有