小波包变换的入门 - - -MATLAB

最近用到小波方面的知识,尤其是小波包变换。

小波包变换的优势:(大部分书上 网上都有,我就简单摘了点过来)
    由于正交小波变换只对信号的低频部分做进一步分解,而对高频部分也即信号的细节部分不再继续分解,所以小波变换能够很好地表征一大类以低频信息为主要成分的信号,但它不能很好地分解和表示包含大量细节信息(细小边缘或纹理)的信号,如非平稳机械振动信号、遥感图象、地震信号和生物医学信号等。与之不同的是,小波包变换可以对高频部分提供更精细的分解,而且这种分解既无冗余,也无疏漏,所以对包含大量中、高频信息的信号能够进行更好的时频局部化分析。

研究了两天,发现如果从头开始研究需要的时间太长,而且如果想真正弄懂小波,还需要了解泛函的知识并且硬着头皮去看那些枯燥的公式。由于我们科研只要明白个大概,能够找到合适的工具来使用,就可以了。因为之前我弄懂傅里叶变换的时候,也是从先会用再到逐渐深入理解的,所以这次我还是先从会用开始。研究了两天之后,发现我小波变换没弄懂什么,小波包先会用了。由于我痛苦的搜了整个网,慢慢理解了一些东西,所以把有用的几个部分拿过来,结合MATLAB,给和我一样想入门的同学一个参考。

首先,小波包的一些基本的基本要弄懂,就是小波包是从原始信号,分级向下分解。如下图所示。
这就是小波包树,其中节点的命名规则是从 (1,0)开始,叫1号, (1,1)是2号,,,,依此类推,(3,0)是7号,(3,7)是14号。 每个节点都有对应的小波包系数,这个系数决定了频率的大小,也就是说频率信息已经有了,但是时域信息在哪里呢? 那就是 order。  这个order就是这些节点的顺序,也就是频率的顺序。
比如,节点的排序是 1,2,3,,,,14, 那么频率就按先1号的频率变化,后2号的,再3号的,,,然后14号的。
小波包变换的入门 <wbr>- <wbr>- <wbr>-MATLAB

                                    图1

来看一个实例:
采样频率为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层,像图1那样,'dmey'使用meyr小波
plot(tt)                 %这个就是画出图1那个图,可以用鼠标在上面点
wpviewcf(tt,1);          %画出时间频率图,如图2

小波包变换的入门 <wbr>- <wbr>- <wbr>-MATLAB

                                   图2

现在开始解释:x轴很简单,就是1024个点,对应1秒,每个点就代表1/1024秒,x轴诡异一下,最后一个数就是1. y轴上显示的数字对应于 图1 中的节点,从下面开始,顺序是 7号节点,8号,10号,9号,,,,11号节点,这个顺序是这么排列的,这是小波包自动排列的,不用管。只要知道怎么查看这个order就可以了。然后,y轴是频率啊,怎么不是 100Hz和300Hz呢? 原因就是MATLAB这里没有显示频率,显示的是order,频率我们要自己算,怎么算呢。我们的采样频率是1024Hz,根据采样定理,奈奎斯特采样频率是512Hz,我们分解了3层,最后一层就是 2^3=8个频率段,每个频率段的频率区间是 512/8=64Hz,对吧,那看图2,颜色重的地方一个是在8那里,一个在13那里,8是第二段,也就是 65-128Hz之间,13是第五段,也就是257-320Hz之间。这样就说通了,正好这个原始信号只有两个频率段,一个100一个300。如果我们不是分解了3层,而是更多层,那么每个频率段包含的频率也就越窄,图上有颜色的地方也会更细,也就是说更精细了,大家可以自己试试。将3改为6试试。由于原始信号的频率在整个1秒钟内都没有改变,所以有颜色的地方是一个横线。

再看一个实例:

有如下的一个信号,该信号的频率从25Hz左右增长到103Hz,信号长度是256,fs就定为256Hz,也是采样1秒。我们用上面的代码来分析这个信号,不过这次分解层数选为4层,也就是有2^4=16个频率段。每个频段是 128/16=8Hz. 这个时候就明白了前面讲到的order的重要性了吧。如果排序不是按照15 16 18 17 21,,23那个顺序拍的,就不可能出现这个随时间而频率增大的图了。从图上还可以看出,频率从第三个,也就是24Hz(3*8),一路走高到13个,也就是13*8=104Hz,正好和信号的图示一样,频率逐渐增大。
小波包变换的入门 <wbr>- <wbr>- <wbr>-MATLAB

小波包变换的入门 <wbr>- <wbr>- <wbr>-MATLAB

小波变换的图像处理%MATLAB2维小波变换经典程序 % FWT_DB.M; % 此示意程序用DWT实现二维小波变换 % 编程时间2004-4-10,编程人沙威 %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc; T=256; % 图像维数 SUB_T=T/2; % 子图维数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 1.调原始图像矩阵 load wbarb; % 下载图像 f=X; % 原始图像 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 2.进行二维小波分解 l=wfilters('db10','l'); % db10(消失矩为10)低通分解滤波器冲击响应(长度为20) L=T-length(l); l_zeros=[l,zeros(1,L)]; % 矩阵行数与输入图像一致,为2的整数幂 h=wfilters('db10','h'); % db10(消失矩为10)高通分解滤波器冲击响应(长度为20) h_zeros=[h,zeros(1,L)]; % 矩阵行数与输入图像一致,为2的整数幂 for i=1:T; % 列变换 row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积FFT row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).'; % 圆周卷积FFT end; for j=1:T; % 行变换 line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ); % 圆周卷积FFT line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ); % 圆周卷积FFT end; decompose_pic=line; % 分解矩阵 % 图像分为四块 lt_pic=decompose_pic(1:SUB_T,1:SUB_T); % 在矩阵左上方为低频分量--fi(x)*fi(y) rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T); % 矩阵右上为--fi(x)*psi(y) lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T); % 矩阵左下为--psi(x)*fi(y) rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T); % 右下方为高频分量--psi(x)*psi(y) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 3.分解结果显示 figure(1); colormap(map); subplot(2,1,1); image(f); % 原始图像 title('original pic'); subplot(2,1,2); image(abs(decompose_pic)); % 分解后图像 title('decomposed pic'); figure(2); colormap(map); subplot(2,2,1); image(abs(lt_pic)); % 左上方为低频分量--fi(x)*fi(y) title('\Phi(x)*\Phi(y)'); subplot(2,2,2); image(abs(rt_pic)); % 矩阵右上为--fi(x)*psi(y) title('\Phi(x)*\Psi(y)'); subplot(2,2,3); image(abs(lb_pic)); % 矩阵左下为--psi(x)*fi(y) title('\Psi(x)*\Phi(y)'); subplot(2,2,4); image(abs(rb_pic)); % 右下方为高频分量--psi(x)*psi(y) title('\Psi(x)*\Psi(y)'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 5.重构源图像及结果显示 % construct_pic=decompose_matrix'*decompose_pic*decompose_matrix; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% l_re=l_zeros(end:-1:1); % 重构低通滤波 l_r=circshift(l_re',1)'; % 位置调整 h_re=h_zeros(end:-1:1); % 重构高通滤波 h_r=circshift(h_re',1)'; % 位置调整 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% top_pic=[lt_pic,rt_pic]; % 图像上半部分 t=0; for i=1:T; % 行插值低频 if (mod(i,2)==0) topll(i,:)=top_pic(t,:); % 偶数行保持 else t=t+1; topll(i,:)=zeros(1,T); % 奇数行为零 end end; for i=1:T; % 列变换 topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圆周卷积FFT end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bottom_pic=[lb_pic,rb_pic]; % 图像下半部分 t=0; for i=1:T; % 行插值高频 if (mod(i,2)==0) bottomlh(i,:)=bottom_pic(t,:); % 偶数行保持 else bottomlh(i,:)=zeros(1,T); % 奇数行为零 t=t+1; end end; for i=1:T; % 列变换 bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圆周卷积FFT end; construct1=bottomch_re+topcl_re; % 列变换重构完毕 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% left_pic=construct1(:,1:SUB_T); % 图像左半部分 t=0; for i=1:T; % 列插值低频 if (mod(i,2)==0) leftll(:,i)=left_pic(:,t); % 偶数列保持 else t=t+1; leftll(:,i)=zeros(T,1); % 奇数列为零 end end; for i=1:T; % 行变换 leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圆周卷积FFT end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% right_pic=construct1(:,SUB_T+1:T); % 图像右半部分 t=0; for i=1:T; % 列插值高频 if (mod(i,2)==0) rightlh(:,i)=right_pic(:,t); % 偶数列保持 else rightlh(:,i)=zeros(T,1); % 奇数列为零 t=t+1; end end; for i=1:T; % 行变换 rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圆周卷积FFT end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% construct_pic=rightch_re+leftcl_re; % 重建全部图像 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 结果显示 figure(3); colormap(map); subplot(2,1,1); image(f); % 源图像显示 title('original pic'); subplot(2,1,2); image(abs(construct_pic)); % 重构源图像显示 title('reconstructed pic'); error=abs(construct_pic-f); % 重构图形与原始图像误值 figure(4); mesh(error); % 误差三维图像 title('absolute error display');
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

怀想天空2011

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值