2代小波示意程序
2維小波變換經典程序
Daubechies小波基的構造
采用多孔trous算法(undecimated wavelet transform)實現小波變換
平移變換平移法(cycle_spinning)消除gibbs效應
提升法97經典程序
消失矩作用的程序
小波插值與小波構造
小波濾波器構造和消噪程序
小波譜分析mallat算法經典程序



2代小波示意程序 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  此程序用提升法實現第二代小波變換
%%  我用的是非整數階小波變換
%%  采用時域實現,步驟先列后行
%%  正變換：分裂，預測，更新；
%%  反變換：更新，預測，合并
%%  只做一層（可以多層，而且每層的預測和更新方程不同）
clear;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  1.調原始圖像矩陣
load wbarb;  %  下載圖像
f=X;         %  原始圖像
% f=[0 0 0 0 0 0 0 0 ;...
%    0 0 0 1 1 0 0 0 ;...
%    0 0 2 4 4 2 0 0 ;...
%    0 1 4 8 8 4 1 0 ;...
%    0 1 4 8 8 4 1 0 ;...
%    0 0 2 4 4 2 0 0 ;...
%    0 0 0 1 1 0 0 0 ;...
%    0 0 0 0 0 0 0 0 ;];  %  原始圖像矩陣
N=length(f);         %  圖像維數
T=N/2;               %  子圖像維數
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%正變換%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   1.列變換
%  A.分裂（奇偶分開）
f1=f([1:2:N-1],:);  %  奇數
f2=f([2:2:N],:);    %  偶數
% f1(:,T+1)=f1(:,1);  %  補列
% f2(T+1,:)=f2(1,:);  %  補行
%  B.預測
for i_hc=1:T;
    high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
end;
% high_frequency_column(T+1,:)=high_frequency_column(1,:);  %  補行
%  C.更新
for i_lc=1:T;
    low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
end;
%  D.合并
f_column([1:1:T],:)=low_frequency_column([1:T],:);
f_column([T+1:1:N],:)=high_frequency_column([1:T],:);
    
    
figure(1)
colormap(map);
image(f);
figure(2)
colormap(map);
image(f_column);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   2.行變換
%  A.分裂（奇偶分開）
f1=f_column(:,[1:2:N-1]);  %  奇數
f2=f_column(:,[2:2:N]);    %  偶數

% f2(:,T+1)=f2(:,1);    %  補行
%  B.預測
for i_hr=1:T;
    high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
end;
% high_frequency_row(:,T+1)=high_frequency_row(:,1);  %  補行
%  C.更新
for i_lr=1:T;
    low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
end;
%  D.合并
f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);
f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);
    
figure(3)
colormap(map);
image(f_row);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%反變換%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   1.行變換
%   A.提取（低頻高頻分開）
f1=f_row(:,[T+1:1:N]);  %  奇數
f2=f_row(:,[1:1:T]);    %  偶數

% f2(:,T+1)=f2(:,1);    %  補行
%  B.更新
for i_lr=1:T;
    low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
end;
%  C.預測
for i_hr=1:T;
    high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
end;
% high_frequency_row(:,T+1)=high_frequency_row(:,1);  %  補行

%  D.合并(奇偶分開合并)
f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);
f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);
    
figure(4)
colormap(map);
image(f_row);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   2.列變換
%  A.提取（低頻高頻分開）
f1=f_row([T+1:1:N],:);  %  奇數
f2=f_row([1:1:T],:);    %  偶數
% f1(:,T+1)=f1(:,1);  %  補列
% f2(T+1,:)=f2(1,:);  %  補行
%  B.更新
for i_lc=1:T;
    low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
end;
%  C.預測
for i_hc=1:T;
    high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
end;
% high_frequency_column(T+1,:)=high_frequency_column(1,:);  %  補行
 
%  D.合并(奇偶分開合并）
f_column([2:2:N],:)=low_frequency_column([1:T],:);
f_column([1:2:N-1],:)=high_frequency_column([1:T],:);
    
    
figure(5)
colormap(map);
image(f_column); 
 


2維小波變換經典程序 
 
%  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');  



Daubechies小波基的構造 
 
%  此程序實現構造小波基
%  periodic_wavelet.m
function ss=periodic_wavelet;
clear;clc;
% global MOMENT;  %  消失矩階數
% global LEFT_SCALET;  %  尺度函數左支撐區(qū)間
% global RIGHT_SCALET;  %  尺度函數右支撐區(qū)間
% global LEFT_BASIS;  %  小波基函數左支撐區(qū)間
% global RIGHT_BASIS;  %  小波基函數右支撐區(qū)間
% global MIN_STEP;  %  最小離散步長
% global LEVEL;  %  計算需要的層數(離散精度）
% global MAX_LEVEL;  %  周期小波最大計算層數

[s2,h]=scale_integer;
[test,h]=scalet_stretch(s2,h);
wave_base=wavelet(test,h);
ss=periodic_waveletbasis(wave_base);
 

function [s2,h]=scale_integer;
%  本函數實現求解小波尺度函數離散整數點的值
%  sacle_integer.m
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區(qū)間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區(qū)間
LEFT_BASIS=1-MOMENT;    %  小波基函數左支撐區(qū)間
RIGHT_BASIS=MOMENT;     %  小波基函數右支撐區(qū)間
MIN_STEP=1/512;          %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度）
MAX_LEVEL=8;  %  周期小波最大計算層數

h=wfilters('db10','r');  %  濾波器系數
h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1;
for i=LEFT_SCALET+1:RIGHT_SCALET-1
    for j=LEFT_SCALET+1:RIGHT_SCALET-1
       k=2*i-j+1;
       if (k>=1&k<=RIGHT_SCALET+1)
       a(i,j)=h(k);  %  矩陣系數矩陣
       else
       a(i,j)=0;
       end
    end
end
[s,w]=eig(a);  %  求特征向量，解的基
s1=s(:,1);
s2=[0;s1/sum(s1);0]; %  根據條件SUM(FI(T))=1,求解;
 
 
%  本函數實現尺度函數經伸縮后的離散值
%  scalet_stretch.m
function [s2,h]=scalet_stretch(s2,h);
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區(qū)間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區(qū)間
LEFT_BASIS=1-MOMENT;  %  小波基函數左支撐區(qū)間
RIGHT_BASIS=MOMENT;  %  小波基函數右支撐區(qū)間
MIN_STEP=1/512;  %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度）
MAX_LEVEL=8;  %  周期小波最大計算層數

for j=1:LEVEL  %  需要計算到尺度函數的層數
   t=0;
   for i=1:2:2*length(s2)-3  %  需要計算的離散點取值（0，1，2，3 -> 1/2, 3/2, 5/2)
      t=t+1;
      fi(t)=0;
      for n=LEFT_SCALET:RIGHT_SCALET;  % 低通濾波器沖擊響應緊支撐判斷
          if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) %  小波尺度函數緊支撐判斷
            fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1);  %  反復應用雙尺度方程求解
          end
      end
   end
   clear s
   n1=length(s2);
   n2=length(fi);
   for i=1:length(s2)+length(fi)  %  變換后的矩陣長度
      if (mod(i,2)==1)
      s(i)=s2((i+1)/2);  %  矩陣奇數下標為小波上一層(0,1,2,3)離散值
      else
      s(i)=fi(i/2);  %  矩陣偶數下標為小波下一層(1/2,3/2,5/2)(經過伸縮變換后)的離散值
      end
   end
   s2=s;
end
 
 

%  采用雙尺度方程求解小波基函數 PSI(T)
%  wavelet.m
function wave_base=wavelet(test,h);
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區(qū)間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區(qū)間
LEFT_BASIS=1-MOMENT;  %  小波基函數左支撐區(qū)間
RIGHT_BASIS=MOMENT;  %  小波基函數右支撐區(qū)間
MIN_STEP=1/512;  %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度）
MAX_LEVEL=8;  %  周期小波最大計算層數

i=0;
for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS;  %  小波基支撐長度 
    s=0;
    for n=1-RIGHT_SCALET:1-LEFT_SCALET  %  g(n)取值范圍
        if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET)  %  尺度函數判斷
          s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1);  %  計算任意精度的小波基函數值  
        end 
    end
    i=i+1;
    wave_base(i)=s;
end
 
采用多孔trous算法(undecimated wavelet transform)實現小波變換 
 
clear;clc;
%%  1.生成信號
f=50;   %  頻率
fs=800; %  采樣率
T=128;  %  信號長度
n=1:T;
y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs));  %  信號
% y=circshift(y.',3).';
%%  2.正變換
l1=wfilters('db4','l')*sqrt(2)/2;     %  參考低通濾波器
l1_zeros=[l1,zeros(1,T-length(l1))];  %  低通濾波器1
h1=wfilters('db4','h')*sqrt(2)/2;     %  參考高通濾波器
h1_zeros=[h1,zeros(1,T-length(h1))];  %  高通濾波器1
low1=ifft(fft(y).*fft(l1_zeros));     %  低頻分量1
high1=ifft(fft(y).*fft(h1_zeros));    %  高頻分量1
l2=dyadup(l1);  %  原濾波器插值
l2_zeros=[l2,zeros(1,T-length(l2))];  %  低通濾波器2
h2=dyadup(h1);  %  原濾波器插值
h2_zeros=[h2,zeros(1,T-length(h2))];  %  高通濾波器2
low2=ifft(fft(low1).*fft(l2_zeros));  %  低頻分量2
high2=ifft(fft(low1).*fft(h2_zeros)); %  高頻分量2
%%  3.反變換
lr2=circshift(l2_zeros(end:-1:1).',1).';  %  重構低通濾波器2
hr2=circshift(h2_zeros(end:-1:1).',1).';  %  重構高通濾波器2
lr1=circshift(l1_zeros(end:-1:1).',1).';  %  重構低通濾波器1
hr1=circshift(h1_zeros(end:-1:1).',1).';  %  重構高通濾波器1
lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2)));  %  重構低頻分量1(lowr=low1)
r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1)));   %  重構源信號(r_s=y)

%%  4.繪圖
figure(1);
plot(y);
title('源信號');
figure(2);
plot(low1,'r');
hold on;
plot(low2,'b');
legend('第一層低頻','第二層低頻');
figure(3);
plot(high1,'r');
hold on;
plot(high2,'b');
legend('第一層高頻','第二層高頻');
figure(4);
plot(low1,'r');
hold on;
plot(lowr,'b.');
legend('第一層低頻','重構第一層低頻');
figure(5);
plot(y,'r');
hold on;
plot(r_s,'b.');
legend('源信號','重構信號');
disp(norm(low1-lowr))
disp(norm(y-r_s)) 



平移變換平移法(cycle_spinning)消除gibbs效應 
 
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  1.原始信號
f=50;    % 信號頻率
fs=800;  % 采樣頻率