小波包變換的原理以及在小目標檢測上的應用
阿新 • • 發佈：2019-01-06
一、利用小波包進行紅外小目標檢測分割：
原理以及程式碼：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！