MATLAB2維小波變換經(jīng)典程序

% FWT_DB.M;
% 此示意程序用DWT實(shí)現(xiàn)二維小波變換
% 編程時(shí)間2004-4-10，編程人沙威
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;
T=256; % 圖像維數(shù)
SUB_T=T/2; % 子圖維數(shù)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.調(diào)原始圖像矩陣
load wbarb; % 下載圖像
f=X; % 原始圖像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.進(jìn)行二維小波分解
l=wfilters('db10','l'); % db10（消失矩為10)低通分解濾波器沖擊響應(yīng)（長(zhǎng)度為20）
L=T-length(l);
l_zeros=[l,zeros(1,L)]; % 矩陣行數(shù)與輸入圖像一致，為2的整數(shù)冪
h=wfilters('db10','h'); % db10（消失矩為10)高通分解濾波器沖擊響應(yīng)（長(zhǎng)度為20）
h_zeros=[h,zeros(1,L)]; % 矩陣行數(shù)與輸入圖像一致，為2的整數(shù)冪
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.分解結(jié)果顯示
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.重構(gòu)源圖像及結(jié)果顯示
% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_re=l_zeros(end:-1:1); % 重構(gòu)低通濾波
l_r=circshift(l_re',1)'; % 位置調(diào)整
h_re=h_zeros(end:-1:1); % 重構(gòu)高通濾波
h_r=circshift(h_re',1)'; % 位置調(diào)整

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
top_pic=[lt_pic,rt_pic]; % 圖像上半部分
t=0;
for i=1:T; % 行插值低頻

if (mod(i,2)==0)
topll(i,:)=top_pic(t,:); % 偶數(shù)行保持
else
t=t+1;
topll(i,:)=zeros(1,T); % 奇數(shù)行為零
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,:); % 偶數(shù)行保持
else
bottomlh(i,:)=zeros(1,T); % 奇數(shù)行為零
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; % 列變換重構(gòu)完畢

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
left_pic=construct1(:,1:SUB_T); % 圖像左半部分
t=0;
for i=1:T; % 列插值低頻

if (mod(i,2)==0)
leftll(:,i)=left_pic(:,t); % 偶數(shù)列保持
else
t=t+1;
leftll(:,i)=zeros(T,1); % 奇數(shù)列為零
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); % 偶數(shù)列保持
else
rightlh(:,i)=zeros(T,1); % 奇數(shù)列為零
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; % 重建全部圖像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 結(jié)果顯示
figure(3);
colormap(map);
subplot(2,1,1);
image(f); % 源圖像顯示
title('original pic');
subplot(2,1,2);
image(abs(construct_pic)); % 重構(gòu)源圖像顯示
title('reconstructed pic');
error=abs(construct_pic-f); % 重構(gòu)圖形與原始圖像誤值
figure(4);
mesh(error); % 誤差三維圖像
title('absolute error display');




clear
clc
%在噪聲環(huán)境下語(yǔ)音信號(hào)的增強(qiáng)
%語(yǔ)音信號(hào)為讀入的聲音文件
%噪聲為正態(tài)隨機(jī)噪聲
sound=wavread('c12345.wav');
count1=length(sound);
noise=0.05*randn(1,count1);
for i=1:count1
signal(i)=sound(i);
end
for i=1:count1
y(i)=signal(i)+noise(i);
end

%在小波基'db3'下進(jìn)行一維離散小波變換
[coefs1,coefs2]=dwt(y,'db3'); %[低頻 高頻]

count2=length(coefs1);
count3=length(coefs2);

energy1=sum((abs(coefs1)).^2);
energy2=sum((abs(coefs2)).^2);
energy3=energy1+energy2;

for i=1:count2
recoefs1(i)=coefs1(i)/energy3;
end
for i=1:count3
recoefs2(i)=coefs2(i)/energy3;
end

%低頻系數(shù)進(jìn)行語(yǔ)音信號(hào)清濁音的判別
zhen=160;
count4=fix(count2/zhen);
for i=1:count4
n=160*(i-1)+1:160+160*(i-1);
s=sound(n);
w=hamming(160);
sw=s.*w;
a=aryule(sw,10);
sw=filter(a,1,sw);
sw=sw/sum(sw);
r=xcorr(sw,'biased');
corr=max(r);
%為清音（unvoice）時(shí)，輸出為1；為濁音（voice）時(shí)，輸出為0
if corr>=0.8
output1(i)=0;
elseif corr<=0.1
output1(i)=1;
end
end
for i=1:count4
n=160*(i-1)+1:160+160*(i-1);
if output1(i)==1
switch abs(recoefs1(i))
case abs(recoefs1(i))<=0.002
recoefs1(i)=0;
case abs(recoefs1(i))>0.002 & abs(recoefs1(i))<=0.003
recoefs1(i)=sgn(recoefs1(i))*(0.003*abs(recoefs1(i))-0.000003)/0.002;
otherwise recoefs1(i)=recoefs1(i);
end
elseif output1(i)==0
recoefs1(i)=recoefs1(i);
end
end

%對(duì)高頻系數(shù)進(jìn)行語(yǔ)音信號(hào)清濁音的判別
count5=fix(count3/zhen);
for i=1:count5
n=160*(i-1)+1:160+160*(i-1);
s=sound(n);
w=hamming(160);
sw=s.*w;
a=aryule(sw,10);
sw=filter(a,1,sw);
sw=sw/sum(sw);
r=xcorr(sw,'biased');
corr=max(r);
%為清音（unvoice）時(shí)，輸出為1；為濁音（voice）時(shí)，輸出為0
if corr>=0.8
output2(i)=0;
elseif corr<=0.1
output2(i)=1;
end
end
for i=1:count5
n=160*(i-1)+1:160+160*(i-1);
if output2(i)==1
switch abs(recoefs2(i))
case abs(recoefs2(i))<=0.002
recoefs2(i)=0;
case abs(recoefs2(i))>0.002 & abs(recoefs2(i))<=0.003
recoefs2(i)=sgn(recoefs2(i))*(0.003*abs(recoefs2(i))-0.000003)/0.002;
otherwise recoefs2(i)=recoefs2(i);
end
elseif output2(i)==0
recoefs2(i)=recoefs2(i);
end
end
%在小波基'db3'下進(jìn)行一維離散小波反變換
output3=idwt(recoefs1, recoefs2,'db3');
%對(duì)輸出信號(hào)抽樣點(diǎn)值進(jìn)行歸一化處理
maxdata=max(output3);
output4=output3/maxdata;
%讀出帶噪語(yǔ)音信號(hào)，存為'101.wav'
wavwrite(y,5500,16,'c101');
%讀出處理后語(yǔ)音信號(hào)，存為'102.wav'
wavwrite(output4,5500,16,'c102');



function [I_W , S] = func_DWT(I, level, Lo_D, Hi_D);
%通過(guò)這個(gè)函數(shù)將I進(jìn)行小波分解，并將分解后的一維向量轉(zhuǎn)換為矩陣形式
% Matlab implementation of SPIHT (without Arithmatic coding stage)
% Wavelet decomposition
% input: I : input image
% level : wavelet decomposition level
% Lo_D : low-pass decomposition filter
% Hi_D : high-pass decomposition filter
% output: I_W : decomposed image vector
% S : corresponding bookkeeping matrix
% please refer wavedec2 function to see more
[C,S] = func_Mywavedec2(I,level,Lo_D,Hi_D);

S(:,3) = S(:,1).*S(:,2); % dim of detail coef nmatrices 求低頻和每個(gè)尺度中高頻的元素個(gè)數(shù)
%st=S(1,3)+S(2,3)*3+S(3,3)*3;%%%%對(duì)前兩層加密
%C(1:st)=0;

L = length(S); %a求S的列數(shù)

I_W = zeros(S(L,1),S(L,2));%設(shè)一個(gè)與原圖像大小相同的全零矩陣

% approx part
I_W( 1:S(1,1) , 1:S(1,2) ) = reshape(C(1:S(1,3)),S(1,1:2)); %將LL層從C中還原為S(1,1)*S(1,2)的矩陣

for k = 2 : L-1 %將C向量中還原出HL,HH,LH 矩陣
rows = [sum(S(1:k-1,1))+1:sum(S(1:k,1))];
columns = [sum(S(1:k-1,2))+1:sum(S(1:k,2))];
% horizontal part
c_start = S(1,3) + 3*sum(S(2:k-1,3)) + 1;
c_stop = S(1,3) + 3*sum(S(2:k-1,3)) + S(k,3);
I_W( 1:S(k,1) , columns ) = reshape( C(c_start:c_stop) , S(k,1:2) );

% vertical part
c_start = S(1,3) + 3*sum(S(2:k-1,3)) + S(k,3) + 1;
c_stop = S(1,3) + 3*sum(S(2:k-1,3)) + 2*S(k,3);
I_W( rows , 1:S(k,2) ) = reshape( C(c_start:c_stop) , S(k,1:2) );

% diagonal part
c_start = S(1,3) + 3*sum(S(2:k-1,3)) + 2*S(k,3) + 1;
c_stop = S(1,3) + 3*sum(S(2:k,3));
I_W( rows , columns ) = reshape( C(c_start:c_stop) , S(k,1:2) );

end

%%%%%%%mallat algorithm%%%%% clc; clear;tic; %%%%original signal%%%% f=100;%%frequence ts=1/800;%%抽樣間隔 N=1:100;%%點(diǎn)數(shù) s=sin(2*ts*pi*f.*N);%%源信號(hào) figure(1) plot(s);%%%源信號(hào)s title('原信號(hào)'); grid on; %%%%小波濾波器%%%% ld=wfilters('db1','l');%%低通 hd=wfilters('db1','h');%%高通 figure(2) stem(ld,'r');%%%低通 grid on; figure(3) stem(hd,'b')%%%高通 grid on; %%%%% tem=conv(s,ld);%%低通和原信號(hào)卷積 ca1=dyaddown(tem);%%抽樣 figure(4) plot(ca1); grid on; tem=conv(s,hd);%%高通和原信號(hào)卷積 cb1=dyaddown(tem);%%抽樣 figure(5) plot(cb1); grid on; %%%%%%%% %[ca3,cb3]=dwt(s,'db1');%%小波變換 %%%%%%%% [lr,hr]=wfilters('db1','r');%%重構(gòu)濾波器 figure(6) stem(lr); figure(7) stem(hr); tem=dyadup(cb1);%%插值 tem=conv(tem,hr);%%卷積 d1=wkeep(tem,100);%%去掉兩頭的分量 %%%%%%%%% tem=dyadup(ca1);%%插值 tem=conv(tem,lr);%%卷積 a1=wkeep(tem,100);%%去掉兩頭的分量 a=a1+d1;%%%重構(gòu)原信號(hào) %%%%%%%%% %a3=idwt(ca3,cb3,'db1',100);%%%小波逆變換 %%%%%%%%% figure(8) plot(a,'.b'); hold on; plot(s,'r'); grid on; title('重構(gòu)信號(hào)和原信號(hào)的比較');toc; %figure(9) %plot(a3,'.b'); %hold on; %plot(s,'r'); %grid on; %title('重構(gòu)信號(hào)和原信號(hào)的比較');
通用函數(shù)

 Allnodes   計(jì)算樹結(jié)點(diǎn)  
appcoef   提取一維小波變換低頻系數(shù)  
appcoef2   提取二維小波分解低頻系數(shù)  
bestlevt   計(jì)算完整最佳小波包樹  
besttree   計(jì)算最佳(優(yōu))樹  
*　 biorfilt   雙正交樣條小波濾波器組  
biorwavf   雙正交樣條小波濾波器  
*　 centfrq   求小波中心頻率  
cgauwavf   Complex Gaussian小波  
cmorwavf   coiflets小波濾波器  
cwt   一維連續(xù)小波變換  
dbaux   Daubechies小波濾波器計(jì)算  
dbwavf   Daubechies小波濾波器   dbwavf(W)    W='dbN'   N=1,2,3,...,50  
ddencmp   獲取默認(rèn)值閾值(軟或硬)熵標(biāo)準(zhǔn)  
depo2ind   將深度-位置結(jié)點(diǎn)形式轉(zhuǎn)化成索引結(jié)點(diǎn)形式  
detcoef   提取一維小波變換高頻系數(shù)  
detcoef2   提取二維小波分解高頻系數(shù)  
disp   顯示文本或矩陣  
drawtree   畫小波包分解樹(GUI)  
dtree   構(gòu)造DTREE類  
dwt   單尺度一維離散小波變換  
dwt2   單尺度二維離散小波變換  
dwtmode   離散小波變換拓展模式  
*　 dyaddown   二元取樣  
*　 dyadup   二元插值  
entrupd   更新小波包的熵值  
fbspwavf   B樣條小波  
gauswavf   Gaussian小波  
get   獲取對(duì)象屬性值  
idwt   單尺度一維離散小波逆變換  
idwt2   單尺度二維離散小波逆變換  
ind2depo   將索引結(jié)點(diǎn)形式轉(zhuǎn)化成深度—位置結(jié)點(diǎn)形式  
*　 intwave   積分小波數(shù)  
isnode   判斷結(jié)點(diǎn)是否存在   
istnode   判斷結(jié)點(diǎn)是否是終結(jié)點(diǎn)并返回排列值  
iswt   一維逆SWT(Stationary Wavelet Transform)變換  
iswt2   二維逆SWT變換  
leaves   　 Determine terminal nodes
mexihat   墨西哥帽小波  meyer   Meyer小波  
meyeraux   Meyer小波輔助函數(shù)  morlet   Morlet小波  
nodease   計(jì)算上溯結(jié)點(diǎn)  
nodedesc   計(jì)算下溯結(jié)點(diǎn)(子結(jié)點(diǎn))  
nodejoin   重組結(jié)點(diǎn)  nodepar   尋找父結(jié)點(diǎn)  
nodesplt   分割(分解)結(jié)點(diǎn)  
noleaves   　 Determine nonterminal nodes
ntnode   　 Number of terminal nodes
ntree   　 Constructor for the class NTREE 
*　 orthfilt   正交小波濾波器組  
plot   繪制向量或矩陣的圖形  
*　 qmf   鏡像二次濾波器  
rbiowavf   　 Reverse biorthogonal spline wavelet filters
read   讀取二進(jìn)制數(shù)據(jù)  readtree   讀取小波包分解樹  
*　 scal2frq   　 Scale to frequency
set   　 
shanwavf   　 Shannon wavelets
swt   一維SWT(Stationary Wavelet Transform)變換  
swt2   二維SWT變換  
symaux   　 Symlet wavelet filter computation.
symwavf   Symlets小波濾波器  
thselect   信號(hào)消噪的閾值選擇  
thodes   　 References
treedpth   求樹的深度  
treeord   求樹結(jié)構(gòu)的叉數(shù)   
upcoef   一維小波分解系數(shù)的直接重構(gòu)  upcoef2   二維小波分解系數(shù)的直接重構(gòu)  
upwlev   單尺度一維小波分解的重構(gòu)  upwlev2   單尺度二維小波分解的重構(gòu)  
wavedec   單尺度一維小波分解  wavedec2   多尺度二維小波分解  
wavedemo   小波工具箱函數(shù)demo  
*　wavefun   小波函數(shù)和尺度函數(shù)  *　 wavefun2   二維小波函數(shù)和尺度函數(shù)  
wavemenu   小波工具箱函數(shù)menu圖形界面調(diào)用函數(shù)  
*　 wavemngr   小波管理函數(shù)  
waverec   多尺度一維小波重構(gòu)  waverec2   多尺度二維小波重構(gòu)  
wbmpen   　 Penalized threshold for wavelet 1-D or 2-D de-noising
wcodemat   對(duì)矩陣進(jìn)行量化編碼  
wdcbm   　 Thresholds for wavelet 1-D using Birge-Massart strategy
wdcbm2   　Thresholds for wavelet 2-D using Birge-Massart strategy 
wden   用小波進(jìn)行一維信號(hào)的消噪或壓縮  
wdencmp   　De-noising or compression using wavelets 
wentropy   計(jì)算小波包的熵  
wextend   　Extend a vector or a matrix 
*　 wfilters   小波濾波器  
wkeep   提取向量或矩陣中的一部分  
*　 wmaxlev   計(jì)算小波分解的最大尺度  
wnoise   產(chǎn)生含噪聲的測(cè)試函數(shù)數(shù)據(jù)  
wnoisest   估計(jì)一維小波的系數(shù)的標(biāo)準(zhǔn)偏差  
wp2wtree   從小波包樹中提取小波樹    　 
wpcoef   計(jì)算小波包系數(shù)  
wpcutree   剪切小波包分解樹  
wpdec   一維小波包的分解  wpdec2   二維小波包的分解  
wpdencmp   用小波包進(jìn)行信號(hào)的消噪或壓縮  
wpfun   小波包函數(shù)  
wpjoin   　重組小波包 
wprcoef   小波包分解系數(shù)的重構(gòu)  
wprec   一維小波包分解的重構(gòu)  wprec2   二維小波包分解的重構(gòu)  
wpsplt   分割（分解）小波包  
wpthcoef   進(jìn)行小波包分解系數(shù)的閾值處理  
wptree   　 顯示小波包樹結(jié)構(gòu)
wpviewcf   　 Plot the colored wavelet packet coefficients. 
wrcoef   對(duì)一維小波系數(shù)進(jìn)行單支重構(gòu)  
wrcoef2   對(duì)二維小波系數(shù)進(jìn)行單支重構(gòu)  
wrev   向量逆序  
write   向緩沖區(qū)內(nèi)存寫進(jìn)數(shù)據(jù)  
wtbo   　 Constructor for the class WTBO 
wthcoef   一維信號(hào)的小波系數(shù)閾值處理  
wthcoef2   二維信號(hào)的小波系數(shù)閾值處理  
wthresh   進(jìn)行軟閾值或硬閾值處理  
wthrmngr   閾值設(shè)置管理  
wtreemgr   管理樹結(jié)構(gòu)