小波包变换的原理以及在小目标检测上的应用

一、利用小波包进行红外小目标检测分割：

原理以及代码：http://blog.csdn.net/chenxieyy/article/details/49180159

http://blog.sina.com.cn/s/blog_8fc890a20101evt5.html

该方法适用于某些图像，通常这类图像在肉眼看起来与周围对比度较大，在目标周围有较强烈的变化，如1 2 6。如果是一些整体比较缓和的图像，则表现不佳，如3 4 。同时附近若有明亮边缘 ，也会造成很大程度的误检，如 8。

function infradDepartion
clear all;
close all;
clc;

tic;
f=imread('E:\A 研究生学习\我的论文~\image\9.jpg');
if ndims(f)>2
    f=rgb2gray(f);  %转为灰度图
end;
imtool(f);
order=2;
depth=4;
f=double(f);
f=medfilt2(f);
T=wpdec2(f,order,'bior4.4'); %2维小波包分解 采用双正交小波bior4.4分解成2层
%plot(T)                 %画出tree图
firstIndex=(order^depth-1)/3;         % 5  depo2ind(order,[depth,0]); %有疑问》公式？还是什么，为什么这样选
lastIndex=((order^(depth+1)-1)/3)-1;  % 9.33 depo2ind(order,[depth,order^depth])-1;
gaus=[];%用来存放节点是不是高斯性的值（是把1存进去，不是把0放进去）.
isGausIndex=[];%用来存放高斯性系数的节点值（也就是 indice value of nodes).
mgausIndex=[];%用来存放gaus数组中等于1的元素下标.
mindex=1;%记录gaus数组中等于1的元素下标
for i=firstIndex:lastIndex
    cfs=wpcoef(T,i);  %求第i个结点的小波包系数 5~9结点
    %cfs=cfs*1.5;   % 改变分解系数矩阵
    igass=judgeGauss(cfs);   %判断是否为高斯系数
    gaus=[gaus,igass];  %由 0 1组成的1维数组
    if igass==1
        isGausIndex=[isGausIndex,i];
        mgausIndex=[mgausIndex,mindex];
    end;
    mindex=mindex+1;
end;
gaus;
pargaus=[];%存放父节点的高斯性
parIndex=[];%存放四个高斯性节点的父节点索引.
%合并具有同一个父节点的四个高斯性节点.
for j=1:order:lastIndex-firstIndex+1
    if j<(lastIndex-firstIndex-2)&gaus(j)==gaus(j+1)&gaus(j+1)==gaus(j+2)&gaus(j+2)==gaus(j+3)&gaus(j+3)==1
        parentIndex=(j+firstIndex-2)/4;
        parIndex=[parIndex,parentIndex];
        cfs=wpcoef(T,parentIndex);
        igass=judgeGauss(cfs);
        if igass==1
            pargaus=[pargaus,igass];
            isGausIndex=[isGausIndex,parentIndex];
        end;
        T=nodejoin(T,parentIndex);
    end;
end;
if numel(parIndex)~=0
    isGausIndex=deleteGausIndex(parIndex,isGausIndex); %这里要把isGausIndex中的高斯性系数的节点值给删除
end;
nGausIndex=numel(isGausIndex);
%抑制低频.
m=firstIndex;
cp=wpcoef(T,m);
T=write(T,'cfs',m,cp);
%将高斯性系数置为0.
for j=1:nGausIndex
    m=isGausIndex(j);
    cp=wpcoef(T,m);
    [h,w]=size(cp);
    y=zeros(h,w);
    T=write(T,'cfs',m,y);
end;

f1=wprec2(T);
means=mean2(f1);
stds=std2(f1);
v=means+stds*3;
[l1,l2]=size(f1);
for i=1:l1
    for j=1:l2
        if f1(i,j)<v
            f1(i,j)=0;
        else
            f1(i,j)=1;
        end;
    end;
end;
figure(1);
subplot(121);
imagesc(f);
title('原始红外图像');
colormap('gray');
subplot(122);
imagesc(uint8(f1));
title('分割后的结果');
toc;

%删除高斯性系数结点
function delGausInd=deleteGausIndex(parentIndex,isGausIndex)
for i=1:numel(parentIndex)
    sonIndex=4*parentIndex(i)+1;
    for j=1:numel(isGausIndex)
        if isGausIndex(j)==sonIndex
            isGausIndex(j:j+3)=[];
            break;
        end;
    end;
    delGausInd=isGausIndex;
end;

%判断小波包系数是否是高斯性系数
function isGauss=judgeGauss(wpacketcoef)
L=numel(wpacketcoef);
[row,col]=size(wpacketcoef);
sum1=0;
sum2=0;
Confidence=0.9;%置信度
for i=1:row
    for j=1:col
        temp1=wpacketcoef(i,j)^4;
        temp2=wpacketcoef(i,j)^2;
        sum1=sum1+temp1;
        sum2=sum2+temp2;
    end;
end;
k=L*(sum1/(sum2^2))-3;
if abs(k)<sqrt(24/(L*(1-Confidence)))
    isGauss=1;
else
    isGauss=0;
end;

以下为程序部分不懂得地方的实验

%小波包基本原理示意图
% close all;
% clear all;
% clc;
% fs=1024;  %采样频率
% f1=100;   %信号的第一个频率
% f2=300;   %信号第二个频率
% t=0:1/fs:1;
% wave=sin(2*pi*f1*t)+sin(2*pi*f2*t);  %生成混合信号
 
 % t=wpdec(wave,3,'dmey');  
% t2 = wpjoin(t,[3;4;5;6]); %将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。
% %因为第一行中小波包树的节点个数是第一层2个，第二层4个，第三层8个。现在将t2就是将第三层的节点再聚合回第二层。
% sNod = read(t,'sizes',[3,4,5,6]); %读取第二层四个节点系数的size
% cfs3  = zeros(sNod(1,:));
% cfs4  = zeros(sNod(2,:));
% cfs5  = zeros(sNod(3,:));
% cfs6  = zeros(sNod(4,:));  %将所有四个节点的小波包系数变为0
% t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);%将四个节点的系数重组到t3小波树中。
% wave2=wprec(t3); %对t3小波树进行重构，获得信号wave2
% figure;
% plot(wave);
% figure;
% plot(wave2);

% %% 时间分辨率
% t=wpdec(wave,1,'dmey');
% figure;
% plot(t);
% wpviewcf(t,1);
% sNod=read(t,'sizes');
% cfs1=zeros(sNod(1,:)); %  注意本行以及以下两行的代码实现~
% cfs1 = cfs1(ones(1,2),:); %长度拓展
% cfs1 = wkeep1(cfs1(:)',1024); %拓展后截取和原来一样的1024个点

二、然后是小波包的具体原理：

小波包分解步骤以及效果图原理：http://blog.sina.com.cn/s/blog_8fc890a20101ecn7.html

综合以上原理自己用matlab画了一个示意图：

以上，每个结点都对应着小波包系数，该系数决定着频率大小，即频域信息。

而时域信息与order有关，整个小波包变换是按照oeder的顺序发生频率变化的！

然后看一下 以上小波包中，顺序排列order的原理：http://blog.sina.com.cn/s/blog_75f0893501015nvb.html

三、小波包分解与重构：http://blog.sina.com.cn/s/blog_8fc890a20101elnd.html

http://www.cnblogs.com/welen/articles/5667217.html（更详细）

主要是代码原理的解释：

例子1：

clear all; 
clc;
fs=1024;  %采样频率
f1=100;   %信号的第一个频率
f2=300;   %信号第二个频率
t=0:1/fs:1;
wave=sin(2*pi*f1*t)+sin(2*pi*f2*t);  %生成混合信号

t=wpdec(wave,3,'dmey'); %小波包分解

% plot(t)                 %画出tree图
% wpviewcf(t,1);     

t2 = wpjoin(t,[3;4;5;6]); %将小波包树的第二行的四个节点收起来，也就是让第二行的节点变为树的最底层节点。
                          %因为第一行中小波包树的节点个数是第一层2个，第二层4个，第三层8个。现在将t2就是将第三层的节点再聚合回第二层。
sNod = read(t,'sizes',[3,4,5,6]); %读取第二层四个节点系数的size
cfs3  = zeros(sNod(1,:));
cfs4  = zeros(sNod(2,:));
cfs5  = zeros(sNod(3,:));
cfs6  = zeros(sNod(4,:));         %将所有四个节点的小波包系数变为0
t3 = write(t2,'cfs',3,cfs3,'cfs',4,cfs4,'cfs',5,cfs5,'cfs',6,cfs6);   %将四个节点的系数重组到t3小波树中。
wave2=wprec(t3);    %对t3小波树进行重构，获得信号wave2
figure;
plot(wave);
figure;
plot(wave2);

上图中可以看出，因为我们把小波树的节点系数都变为0了，所以信号也就全为0了，所以wave2是一个0向量。

进一步，如果我们只聚合第二层中的某几个节点，比如 4和5，那么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 用 meyr小波进行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);

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

OK！