小波变换和小波阈值法去噪
 
1. 小波变换
 
小波变换是一种信号的时间——尺度（时间——频率）分析方法，它具有多分辨分析的特点，而且在时频两域都具有表征信号局部特征的能力，是一种窗口大小固定不变但其形状可改变，时间窗和频率窗都可以改变的时频局部化分析方法。即在低频部分具有较低的时间分辨率和较高的频率分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于分析非平稳的信号和提取信号的局部特征，所以小波变换被誉为分析处理信号的显微镜。
 
傅里叶是将信号分解成一系列不同频率的正余弦函数的叠加，同样小波变换是将信号分解为一系列的小波函数的叠加(或者说不同尺度、时间的小波函数拟合)，而这些小波函数都是一个母小波经过平移和尺度伸缩得来的，如下图。
 
 
小波变换常见的形式有连续小波变换(CWT)、离散小波变换(DWT)等。连续小波变换是在尺度基础上连续变换的，做信号的小波分析得到的是幅值，a时间的三维图，对应的a值所截得的曲线即为该尺度的小波图形。而离散小波变换常用的是二进小波变换，对尺度和时间进行离散化处理。
 
 
 
CWT连续小波变换
 
 
 
 
CWT步骤：
 
首先选择一个小波基函数，固定一个尺度因子，将它与信号的初始段进行比较；
 通过CWT的计算公式计算小波系数（反映了当前尺度下的小波与所对应的信号段的相似程度）；
 改变平移因子，使小波沿时间轴位移，重复上述两个步骤完成一次分析；
 增加尺度因子，重复上述三个步骤进行第二次分析；
 循环执行上述四个步骤，直到满足分析要求为止。
 
 
 
 
 

 
 
 
 
 连续小波变换是在尺度基础上连续变换的，做信号的小波分析得到的是幅值，a时间的三维图，对应的a值所截得的曲线即为该尺度的小波图形。而集散小波变换常用的是二进小波变换。
 
 
    但是，cwt的结果都相当于DWT中的细节信息（即所谓DWT中的高频信息。虽然越向后频率越低，有时已不能用“高频”来形容了，但这时的高频是相对概念，是相对于同阶逼近信息还是高的），只是其尺度是连续的尺度越大频率越低，一直低下去。
 
 
    morlet等小波只能做CWT，有些是因为没法儿构造尺度函数，有些是根本就没有逆变换（只有满足某些条件，CWT才存在逆变换，这与小波基有关），有些是如何离散化也不能构成正交或双正交基，甚至按照二进制的离散化不能构成紧支的框架，所以它们通常不能做DWT，也就没有逆变换、重构一说了。
 DWT离散小波变换
 

 离散小波变换DWT对尺度参数按幂级数进行离散化处理，对时间进行均匀离散化取值如二进制离散化尺度时间为2,4,6,8...2n（要求采样率满足尼奎斯特采样定理），常用于信号的多分辨分析、信号分解重构。
 

 
 
 
 

 多分辨分析也称为多尺度分析，是建立在函数空间概念上的理论。在不同的尺度和时间下，分别构造了尺度函数向量组合小波函数向量组，也即是尺度函数向量空间V与小波函数向量空间W，在一定层次下，信号在尺度空间做卷积所得到的是信号的近似、低频信息，信号在小波空间W做卷积所得到的是信号的细节、高频信息。(注意：尺度与分解层数不是一个概念，尺度与频率成反比的，分解层数是对频率的范围进行一定的划分)。
 
 
 
 
 
 
 

  在多分辨分析中，如正交小波变换可以等效为一组镜像滤波的过程，即信号通过一个分解高通滤波器和分解低通滤波器，自然的高通滤波器输出对应的信号的高频分量部分，称为细节分量，低通滤波器输出对应了信号的相对较低的频率分量部分，称为近似分量。对应的快速算法称为Mallat算法。
 
   
 
 
 
 

 
 

  
 

 小波分解重构过程（其中CA为低频信息、近似分量，CD为高频、细节分量）：
 

  
 

 
 
 
 
小波阈值去噪
 
 
通常情况下， 我们在从设备上采集到的信号都是具有一定的噪声的，大多数情况下，可认为这种噪声为高斯白噪声。被噪声污染的信号=干净的信号+噪声。
  为什么要使用阈值：由于信号在空间上(或者时间域)是有一定连续性的，因此在小波域，有效信号所产生的小波系数其模值往往较大；而高斯白噪声在空间上(或者时间域)是没有连续性的，因此噪声经过小波变换，在小波阈仍然表现为很强的随机性，通常仍认为是高斯白噪的。 那么就得到这样一个结论：在小波域，有效信号对应的系数很大，而噪声对应的系数很小。 刚刚已经说了，噪声在小波域对应的系数仍满足高斯白噪分布。如果在小波域，噪声的小波系数对应的方差为sigma，那么根据高斯分布的特性，绝大部分(99.99%)噪声系数都位于[-3*sigma，3*sigma]区间内(切比雪夫不等式， 3sigma准则)。因此，只要将区间[-3*sigma，3*sigma]内的系数置零(这就是常用的硬阈值函数的作用)，就能最大程度抑制噪声的，同时只是稍微损伤有效信号。将经过阈值处理后的小波系数重构，就可以得到去噪后的信号。 常用的软阈值函数，是为了解决硬阈值函数“一刀切”导致的影响(模小于3*sigma的小波系数全部切除，大于3*sigma全部保留，势必会在小波域产生突变，导致去噪后结果产生局部的抖动，类似于傅立叶变换中频域的阶跃会在时域产生拖尾)。软阈值函数将模小于3*sigma的小波系数全部置零，而将模大于3*sigma的做一个比较特殊的处理，大于3*sigma的小波系数统一减去3*sigma，小于-3*sigma的小波系数统一加3*sigma。经过软阈值函数的作用，小波系数在小波域就比较光滑了，因此用软阈值去噪得到的图象看起来很平滑，类似于冬天通过窗户看外面一样，像有层雾罩在图像上似的。 
 比较硬阈值函数去噪和软阈值函数去噪：硬阈值函数去噪所得到的峰值信噪比(PSNR)较高，但是有局部抖动的现象；软阈值函数去噪所得到的PSNR不如硬阈值函数去噪，但是结果看起来很平滑，原因就是软阈值函数对小波系数进行了较大的 “社会主义改造”，小波系数改变很大。因此各种各样的阈值函数就出现了，其目的我认为就是要使大的系数保留，小的系数被剔出，而且在小波域系数过渡要平滑。
 如何估计小波域噪声方差sigma的估计，这个很简单：把信号做小波变换，在每一个子带利用robust estimator估计就可以(可能高频带和低频带的方差不同)。 robust estimator就是将子带内的小波系数模按大小排列，然后取最中间那个，然后把最中间这个除以0.6745就得到噪声在某个子带内的方差sigma。利用这个sigma，然后选种阈值函数，就可以去去噪了，在matlab有实现api可使用。
 
 
小波阈值去噪过程
 
 

  
 

 
 

 在小波分析中经常用到近似和细节，近似表示信号的高尺度，即低频信息；细节表示信号的低尺度，即高频信息。对含有噪声的信号，噪声分量的主要能量集中在小波解的细节分量中。
 
 在以上过程中，小波基和分解层数的选择，阈值的选取规则，和阈值函数的设计，都是影响最终去噪效果的关键因素。
 

  
 

 1、小波基的选择
 

  
 

 可参考 http://blog.csdn.net/jbb0523/article/details/42586749 博文，一般选取小波基函数要从支撑长度、消失矩、对称性、正则性以及相似性等进行综合考虑。由于小波基函数在处理信号时各有特点，且没有任何一种小波基函数可以对所有类型信号都取得最优的去噪效果。一般来讲，db小波系和sym小波系在语音去噪中是经常会被用到的两族小波基。
 

 2、分解层数的选择
 

 对于一个要采集的信号，根据奈奎斯采样定理，其采样频率>= 2*信号的最大频率。而其他噪声频率如高斯白噪声的信号是幅度分布服从高斯分布，功率谱密度服从均匀分布的，并且与有效信号进行混合叠加的。
 

 在小波分解中，分解层数的选择也是非常重要的一步。取得越大，则噪声和信号表现的不同特性越明显，越有利于二者的分离。但另一方面，分解层数越大，重构到的信号失真也会越大，在一定程度上又会影响最终去噪的效果。因此在应用时要格外注意处理好两者之间的矛盾，选择一个合适的分解尺度。
 
 通常小波分解的频段范围与采样频率有关。若N层分解，则各个频段大小为Fs/2/2^N 。例如：一个原始信号，经历的时间长度为2秒，采样了2000个点，那么做除法，可得出采样频率为1000hz，由采样定理（做除法）得该信号的最大频率为500hz，那么对该信号做3层的DWT，一阶细节的频段为250－500hz，一阶逼近的频段为小于250hz，二阶细节的频段为125－250hz，逼近的频段为小于125hz，三阶细节的频段约为62.5－125hz，逼近的频段为小于62.5hz。对于更多阶的分解也是以此类推的。
 

 3、阈值的选取
 

  
 

 在小波域，有效信号对应的系数很大，而噪声对应的系数很小。噪声在小波域对应的系数仍满足高斯白噪分布。
 

 阈值选择规则基于模型 y = f(t) + e，e是高斯白噪声N(0,1)。因此可以通过小波系数、或者原始信号来进行评估能够消除噪声在小波域的阈值。
 

 目前常见的阈值选择方法有：固定阈值估计、极值阈值估计、无偏似然估计以及启发式估计等（N为信号长度）。
 

 
 

 
 

 
 

 
 

  
 

  
 

 
 

  
 

 一般来讲，极值阈值估计和无偏似然估计方法比较保守，当噪声在信号的高频段分布较少时，这两种阈值估计方法效果较好可以将微弱的信号提取出来。而固定阈值估计和启发式阈值估计去噪比较彻底，在去噪时显得更为有效，但是也容易把有用的信号误认为噪声去掉。
 

  
 

 4、 阈值函数选择
 

  
 

 确定了高斯白噪声在小波系数(域)的阈值门限之后，就需要有个阈值函数对这个含有噪声系数的小波系数进行过滤，去除高斯噪声系数，常用的阈值函数有软阈值和硬阈值方法，很多文献论文中也有在阈值函数进行一些大量的改进和优化。
 

 
 

 
 

  
 

 软硬阈值函数优缺点对比：
 

  
 

  
 

  
 

 硬阈值函数在均方误差意义上优于软阈值法，但是信号会产生附加震荡，产生跳跃点，不具有原始信号的平滑性。
 

 软阈值估计得到的小波系数整体连续性较好，从而使估计信号不会产生附加震荡，但是优于会压缩信号，会产生一定的偏差，
 直接影响到重构的信号与真实信号的逼近程度。
 

  
 

 5、 matlab中小波工具箱
 

  
 

 
 

 
 

  

摘自：http://zh.wikipedia.org/wiki/%E5%B0%8F%E6%B3%A2%E5%88%86%E6%9E%90
 
定义：wavelet是指小型波(在傅里叶分析里的弦波是大型波)，简单说来，小波(wavelet)是一个衰减迅速的振荡。
 
小波函数：小波只有时域表示，作为小波函数。例如墨西哥帽小波。
 
 
     
 
 
小波历史：
 1909年: 第一个小波（Haar小波）由Alfred Haar提出
 1981年: Jean Morlet提出小波概念
 1984年: Morlet和 Grossman发明了"wavelet"一词
 1985年: Meyer小波
 1987年: 第一届国际研讨会在法国
 1988年: Stephane Mallat和 Meyer 提出了多清晰度的概念
 1988年: 多贝西小波
 1989年: 快速小波转换 （FWT）
 1990年代: 离散小波变换 （DWT）
 1999年: 方向小波转换
 2000年: JPEG2000
 
小波转换的优点：
 
 
1. 可以同时观察频率和时间轴，在频率高时有较好的时间清晰度，在频率低时有较好的频率清晰度。
 2. 有快速小波转换可以加速运算。
 3. 可以分离出信号的精细或粗糙成分。
 4. 在小波理论中，可以用较少的小波系数去逼近一个函数。
 5. 对信号去噪或压缩信号时，不会对信号造成明显的破坏。
 
 
小波变换的实质: 
 
小波变换的公式有内积形式和卷积形式，两种形式的实质都是一样的．它要求的就是一个个小波分量的系数也就是＂权＂．其直观意义就是首先用一个时窗最窄，频窗最宽的小波作为尺子去一步步地＂量＂信号，也就是去比较信号与小波的相似程度．信号局部与小波越相似，则小波变换的值越大，否则越小！当一步比较完成后，再将尺子拉长一倍，又去一步步地比较，从而得出一组组数据．如此这般循环，最后得出的就是信号的小波分解（小波级数）．当然这只是一种粗略的解释．
 

 
小波分析、小波降噪matlab代码实现
 
软阈值、硬阈值、固定阈值三种方式
1. 简介
2. 操作步骤
3. 直接上代码
4. 运行结果显示
1. 噪声信号图像
2. 硬阈值去噪图像
3. 软阈值去噪图像
4. 固定阈值去噪图像
5. 去噪结果的误差对比
   
  
 

 
软阈值、硬阈值、固定阈值三种方式
 
1. 简介
 
所谓的小波就是指微小的具有波动性的波形，从信号学角度来看，是一个信号滤波的过程。目前，小波去噪已然是一种喜闻乐见的噪声处理方式。
 本文使用matlab进行小波降噪，采用软阈值、硬阈值、和固定阈值三种方式，对一维噪声数据进行降噪处理，并输出误差SNR,RMSE,以及降噪前后的噪声图像 。
 
2. 操作步骤
 
主要有三个步骤：
 
对被噪声污染的信号进行小波变换处理；
通过对变换后的小波系数进行处理，去除其中的噪声
小波逆变换，得到去噪后的信号
 
3. 直接上代码
 
%%初始化程序
clear,clc
t1=clock;
 %% 载入噪声信号数据，数据为.mat格式，并且和程序放置在同一个文件夹下
load('filename.mat');%matrix
YSJ= filename;
 %% 数据预处理，数据可能是存储在矩阵或者是EXCEL中的二维数据，衔接为一维的，如果数据是一维数据，此步骤也不会影响数据
[c,l]=size(YSJ);
Y=[];
for i=1:c
    Y=[Y,YSJ(i,:)];
end
[c1,l1]=size(Y);
X=[1:l1];
 %% 绘制噪声信号图像
figure(1);
plot(X,Y);
xlabel('横坐标');
ylabel('纵坐标');
title('原始信号');
 %% 硬阈值处理
lev=3;
xd=wden(Y,'heursure','h','one',lev,'db4');%硬阈值去噪处理后的信号序列
figure(2)
plot(X,xd)
xlabel('横坐标');
ylabel('纵坐标');
title('硬阈值去噪处理')
set(gcf,'Color',[1 1 1])
 %% 软阈值处理
lev=3;
xs=wden(Y,'heursure','s','one',lev,'db4');%软阈值去噪处理后的信号序列
figure(3)
plot(X,xs)
xlabel('横坐标');
ylabel('纵坐标');
title('软阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 固定阈值后的去噪处理
lev=3;
xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定阈值去噪处理后的信号序列
figure(4)
plot(X,xz);
xlabel('横坐标');
ylabel('纵坐标');
title('固定阈值后的去噪处理')
set(gcf,'Color',[1 1 1])
%% 计算信噪比SNR
Psig=sum(Y*Y')/l1;
Pnoi1=sum((Y-xd)*(Y-xd)')/l1;
Pnoi2=sum((Y-xs)*(Y-xs)')/l1;
Pnoi3=sum((Y-xz)*(Y-xz)')/l1;
SNR1=10*log10(Psig/Pnoi1);
SNR2=10*log10(Psig/Pnoi2);
SNR3=10*log10(Psig/Pnoi3);
%% 计算均方根误差RMSE
RMSE1=sqrt(Pnoi1);
RMSE2=sqrt(Pnoi2);
RMSE3=sqrt(Pnoi3);
%% 输出结果
disp('-------------三种阈值设定方式的降噪处理结果---------------'); 
disp(['硬阈值去噪处理的SNR=',num2str(SNR1),'，RMSE=',num2str(RMSE1)]);
disp(['软阈值去噪处理的SNR=',num2str(SNR2),'，RMSE=',num2str(RMSE2)]);
disp(['固定阈值后的去噪处理SNR=',num2str(SNR3),'，RMSE=',num2str(RMSE3)]);
t2=clock;
tim=etime(t2,t1);
disp(['------------------运行耗时',num2str(tim),'秒-------------------'])
 
4. 运行结果显示
 
1. 噪声信号图像
 
 
2. 硬阈值去噪图像
 
 
3. 软阈值去噪图像
 
 
4. 固定阈值去噪图像
 
 
5. 去噪结果的误差对比
 

 由于本人能力有限，欢迎大家指出不足之处，我将尽力完善代码

        小波与小波包、小波包分解与信号重构、小波包能量特征提取   (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