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');
clear
clc
%在噪聲環(huán)境下語音信號的增強
%語音信號為讀入的聲音文件
%噪聲為正態(tài)隨機噪聲
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'下進行一維離散小波變換
[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
%低頻系數進行語音信號清濁音的判別
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）時，輸出為1；為濁音（voice）時，輸出為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
%對高頻系數進行語音信號清濁音的判別
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）時，輸出為1；為濁音（voice）時，輸出為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'下進行一維離散小波反變換
output3=idwt(recoefs1, recoefs2,'db3');
%對輸出信號抽樣點值進行歸一化處理
maxdata=max(output3);
output4=output3/maxdata;
%讀出帶噪語音信號，存為'101.wav'
wavwrite(y,5500,16,'c101');
%讀出處理后語音信號，存為'102.wav'
wavwrite(output4,5500,16,'c102');
function [I_W , S] = func_DWT(I, level, Lo_D, Hi_D);
%通過這個函數將I進行小波分解，并將分解后的一維向量轉換為矩陣形式
% 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 求低頻和每個尺度中高頻的元素個數
%st=S(1,3)+S(2,3)*3+S(3,3)*3;%%%%對前兩層加密
%C(1:st)=0;
L = length(S); %a求S的列數
I_W = zeros(S(L,1),S(L,2));%設一個與原圖像大小相同的全零矩陣
% 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;%%點數 s=sin(2*ts*pi*f.*N);%%源信號 figure(1) plot(s);%%%源信號s title('原信號'); 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);%%低通和原信號卷積 ca1=dyaddown(tem);%%抽樣 figure(4) plot(ca1); grid on; tem=conv(s,hd);%%高通和原信號卷積 cb1=dyaddown(tem);%%抽樣 figure(5) plot(cb1); grid on; %%%%%%%% %[ca3,cb3]=dwt(s,'db1');%%小波變換 %%%%%%%% [lr,hr]=wfilters('db1','r');%%重構濾波器 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;%%%重構原信號 %%%%%%%%% %a3=idwt(ca3,cb3,'db1',100);%%%小波逆變換 %%%%%%%%% figure(8) plot(a,'.b'); hold on; plot(s,'r'); grid on; title('重構信號和原信號的比較');toc; %figure(9) %plot(a3,'.b'); %hold on; %plot(s,'r'); %grid on; %title('重構信號和原信號的比較');
通用函數
Allnodes 計算樹結點 appcoef 提取一維小波變換低頻系數 appcoef2 提取二維小波分解低頻系數 bestlevt 計算完整最佳小波包樹 besttree 計算最佳(優(yōu))樹 * biorfilt 雙正交樣條小波濾波器組 biorwavf 雙正交樣條小波濾波器 * centfrq 求小波中心頻率 cgauwavf Complex Gaussian小波 cmorwavf coiflets小波濾波器 cwt 一維連續(xù)小波變換 dbaux Daubechies小波濾波器計算 dbwavf Daubechies小波濾波器 dbwavf(W) W='dbN' N=1,2,3,...,50 ddencmp 獲取默認值閾值(軟或硬)熵標準 depo2ind 將深度-位置結點形式轉化成索引結點形式 detcoef 提取一維小波變換高頻系數 detcoef2 提取二維小波分解高頻系數 disp 顯示文本或矩陣 drawtree 畫小波包分解樹(GUI) dtree 構造DTREE類 dwt 單尺度一維離散小波變換 dwt2 單尺度二維離散小波變換 dwtmode 離散小波變換拓展模式 * dyaddown 二元取樣 * dyadup 二元插值 entrupd 更新小波包的熵值 fbspwavf B樣條小波 gauswavf Gaussian小波 get 獲取對象屬性值 idwt 單尺度一維離散小波逆變換 idwt2 單尺度二維離散小波逆變換 ind2depo 將索引結點形式轉化成深度—位置結點形式 * intwave 積分小波數 isnode 判斷結點是否存在 istnode 判斷結點是否是終結點并返回排列值 iswt 一維逆SWT(Stationary Wavelet Transform)變換 iswt2 二維逆SWT變換 leaves Determine terminal nodes mexihat 墨西哥帽小波 meyer Meyer小波 meyeraux Meyer小波輔助函數 morlet Morlet小波 nodease 計算上溯結點 nodedesc 計算下溯結點(子結點) nodejoin 重組結點 nodepar 尋找父結點 nodesplt 分割(分解)結點 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 讀取二進制數據 readtree 讀取小波包分解樹 * scal2frq Scale to frequency set shanwavf Shannon wavelets swt 一維SWT(Stationary Wavelet Transform)變換 swt2 二維SWT變換 symaux Symlet wavelet filter computation. symwavf Symlets小波濾波器 thselect 信號消噪的閾值選擇 thodes References treedpth 求樹的深度 treeord 求樹結構的叉數 upcoef 一維小波分解系數的直接重構 upcoef2 二維小波分解系數的直接重構 upwlev 單尺度一維小波分解的重構 upwlev2 單尺度二維小波分解的重構 wavedec 單尺度一維小波分解 wavedec2 多尺度二維小波分解 wavedemo 小波工具箱函數demo * wavefun 小波函數和尺度函數 * wavefun2 二維小波函數和尺度函數 wavemenu 小波工具箱函數menu圖形界面調用函數 * wavemngr 小波管理函數 waverec 多尺度一維小波重構 waverec2 多尺度二維小波重構 wbmpen Penalized threshold for wavelet 1-D or 2-D de-noising wcodemat 對矩陣進行量化編碼 wdcbm Thresholds for wavelet 1-D using Birge-Massart strategy wdcbm2 Thresholds for wavelet 2-D using Birge-Massart strategy wden 用小波進行一維信號的消噪或壓縮 wdencmp De-noising or compression using wavelets wentropy 計算小波包的熵 wextend Extend a vector or a matrix * wfilters 小波濾波器 wkeep 提取向量或矩陣中的一部分 * wmaxlev 計算小波分解的最大尺度 wnoise 產生含噪聲的測試函數數據 wnoisest 估計一維小波的系數的標準偏差 wp2wtree 從小波包樹中提取小波樹 wpcoef 計算小波包系數 wpcutree 剪切小波包分解樹 wpdec 一維小波包的分解 wpdec2 二維小波包的分解 wpdencmp 用小波包進行信號的消噪或壓縮 wpfun 小波包函數 wpjoin 重組小波包 wprcoef 小波包分解系數的重構 wprec 一維小波包分解的重構 wprec2 二維小波包分解的重構 wpsplt 分割（分解）小波包 wpthcoef 進行小波包分解系數的閾值處理 wptree 顯示小波包樹結構 wpviewcf Plot the colored wavelet packet coefficients. wrcoef 對一維小波系數進行單支重構 wrcoef2 對二維小波系數進行單支重構 wrev 向量逆序 write 向緩沖區(qū)內存寫進數據 wtbo Constructor for the class WTBO wthcoef 一維信號的小波系數閾值處理 wthcoef2 二維信號的小波系數閾值處理 wthresh 進行軟閾值或硬閾值處理 wthrmngr 閾值設置管理 wtreemgr 管理樹結構
您的評論需要經過審核才能顯示
有用
有用
有用