N=128;   % 采樣點 
%  信號賦值
n=1:N;
y=sin(2*pi*f*n/fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  2.噪聲
noise=0.4*rand(1,128);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  3.染噪信號
y_noise=y+noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  4.硬消噪采用cycle_spinning技術
%  累加量
z5=zeros(1,N);
%  平移變換頻移法
for i=1:N;
    z=circshift(y_noise.',i-1).';  %  源信號右平移
    [z1,z2]=lwt(z,'db3');    %  小波正變換
    z2=zeros(1,N/2);         %  高頻分量全部為零（主要噪聲，硬消噪）
    z3=ilwt(z1,z2,'db3');    %  小波反變換  
    z4=circshift(z3.',-(i-1)).';  %  變換后信號左平移  
    z5=z5+z4/N;  %  平均  
end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  5.顯示
error=norm(y-z5)/norm(y);  % 相對誤差
figure(1);
subplot(2,1,1)
plot(y,'r');
legend('源信號');
subplot(2,1,2);
plot(y_noise);
legend('染噪信號');

figure(2);
subplot(2,1,1)
plot(y,'r');
legend('源信號');
title(error);
subplot(2,1,2);
plot(z5);
legend('消噪后信號'); 



提升法97經(jīng)典程序 
 
%%  本程序?qū)崿F(xiàn)任意偶數(shù)大小圖像第二代雙正交97提升小波變換
%%  注1： 采用標準正交方法，對行列采用不同矩陣（和matlab里不同）
%%  注2： 為了保證正交，所有邊界處理，全部采用循環(huán)處理
%%  注3： 正交性驗證，將單位陣帶入函數(shù)，輸出仍是單位陣（matlab不具有此性質(zhì)）
%%  注4： 此程序是矩陣實現(xiàn)，所以圖像水平分量和垂直分量估計被交換位置
%%  注5： 此程序?qū)崿F(xiàn)的是類小波（wavelet-like）變換，是介于小波包變換與小波變換之間的變換
%%  注6： 此程序每層變換相對原圖像矩陣，產(chǎn)生的矩陣都是正交陣，這和小波包一致
%%  注7： 但小波變換每層產(chǎn)生的矩陣，是相對每個待分解子塊的正交矩陣，而不是原圖像的正交矩陣
%%  注8： 且小波變換產(chǎn)生的正交矩陣維數(shù)，隨分解層數(shù)2分減少
%%  注9:  提升系數(shù)可以在MATLAB7.0以上版本，用liftwave('9.7')獲取，這里直接給出，考慮兼容性
%%  注10：由于MATLAB數(shù)組下標從1開始，所以注意奇偶序列的變化
%%  注11：d為對偶上升，即預測；p為原上升，即更新
%%  編程人  沙威 安徽大學
%%  編程時間  2004年12月18日
%%  x輸入圖像,y輸出圖像
%%  flag_trans為正變換或反變換標志，0執(zhí)行正變換，1執(zhí)行反變換
%%  flag_max,是否最大層數(shù)變換標志,0執(zhí)行用戶設定層數(shù),1執(zhí)行最大層數(shù)變換
%%  layer,用戶層數(shù)設置（小于最大層）
function y=db97(x,flag_trans,flag_max,layer);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  1.輸入?yún)?shù)檢查
%  矩陣維數(shù)判斷
[sa,sb]=size(x);
if (sa~=sb)       %  防止非圖像數(shù)據(jù)
    errordlg('非圖像數(shù)據(jù)！');
    error('非圖像數(shù)據(jù)！');
end;
%  變換標志判斷
[sa,sb]=size(flag_trans);
if ((sa~=1) | (sb~=1))                  %  變換標志錯誤
    errordlg('變換標志錯誤！');
    error('變換標志錯誤！');
end;
if ((flag_trans~=1) & (flag_trans~=0))  %  變換標志錯誤
    errordlg('變換標志錯誤！');
    error('變換標志錯誤！');
end;
%  最大層數(shù)標志判斷
[sa,sb]=size(flag_max);
if ((sa~=1) | (sb~=1))              %  最大層數(shù)標志錯誤
    errordlg('最大層數(shù)標志錯誤！');
    error('最大層數(shù)標志錯誤！');
end;
if ((flag_max~=1) & (flag_max~=0))  %  最大層數(shù)標志錯誤
    errordlg('最大層數(shù)標志錯誤！');
    error('最大層數(shù)標志錯誤！');
end;
%  用戶設置層數(shù)判斷
if (flag_max~=1)
    [sa,sb]=size(layer);
    if ((sa~=1) | (sb~=1))       %  層數(shù)設置錯誤
        errordlg('層數(shù)設置錯誤！');
        error('層數(shù)設置錯誤！');
    end;
    if (flag_max<0)              %  層數(shù)設置錯誤
        errordlg('層數(shù)設置錯誤！');
        error('層數(shù)設置錯誤！');
    end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  2.提升系數(shù)確定
%  t1=liftwave('9.7');  %  獲取提升系數(shù)(MATLAB7.0以后)
d1=[-1.586100000000000e+000,-1.586134342069360e+000];
p1=[1.079600000000000e+000,-5.298011857188560e-002];
d2=[-8.829110755411875e-001,-8.829110755411875e-001];
p2=[4.435068520511142e-001,1.576123746148364e+000];
d3=-8.698644516247808e-001;
p3=-1.149604398860242e+000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  3.分解層數(shù)確定
%  采用用戶輸入和自動給出最大層數(shù)兩種方法
N=length(x); %  矩陣大小
S=N;         %  變量
s=log2(N);   %  最大循環(huán)次數(shù)
n1=N/2;      %  初始一半矩陣大小
n2=N;        %  初始矩陣大小
u=0;         %  初始值
%  對非2的整數(shù)冪大小圖像確定最大分解層數(shù)
for ss=1:s
    if (mod(S,2)==0)
        u=u+1;
        S=S/2;
    end;
end;
u=u-1;            %  分解最大層數(shù)減1(后面的邊界處理造成)
%  最大層數(shù)確定
if (flag_max==0)  %  手動輸入
    T=layer;      %  用戶輸入值
else              %  自動確定最大層數(shù)
    T=u;          %  分解最大層數(shù)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  4.最大層數(shù)和圖像大小檢查
if (T>u)          %  防止用戶層數(shù)越界
    errordlg('已超過最大分解層數(shù)！或者非偶數(shù)大小圖像！');
    error('已超過最大分解層數(shù)！或者非偶數(shù)大小圖像！');
end;
if (mod(N,2)~=0)  %  防止圖像大小錯誤
    errordlg('非偶數(shù)大小圖像！');
    error('非偶數(shù)大小圖像！');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  5.提升法正變換
if (flag_trans==0)
    for time=1:T;
        %  行正變換
        
        % d;
        x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);  
        x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:); 
        
        % p;
        x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:); 
        x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)*x1([2:n1],:); 
        x([n1+1:n2],:)=x1([1:n1],:);
 
        % d;        
        x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:);
        x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1)*x([2:n1],:);
        
        % p;
        x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:);
        x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:);
        x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1)*x([n1+3:n2],:);
        
        % 歸一
        x([1:n1],:)=p3*x([1:n1],:);
        x([n1+1:n2],:)=d3*x([n1+1:n2],:);
        clear x1;
        
        %  列正變換
        
        % d;
        x1(:,[1:n1])=x(:,[2:2:n2]);
        
        % p;
        x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);  
        x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1(:,[2:n1]); 
        x(:,[n1+1:n2])=x1(:,[1:n1]);
        
        % d;
        x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);
        x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p1(2)*x(:,[2:n1]);
 
        % p;        
        x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);
        x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2)*x(:,[n1+2:n2]);
        
        % d;
        x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1);
        x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);
        x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2)*x(:,[2:n1-1]);
        
        % 歸一
        x(:,[1:n1])=d3*x(:,[1:n1]);
        x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]);
        clear x1;
        
        n2=n2/2; %  原大小
        n1=n2/2; %  一半大小
     end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%  6.提升法反變換
else
    n2=N/(2.^(T-1)); %  分解最小子塊維數(shù)
    n1=n2/2;
    for time=1:T;
        %  行反變換
        
        %  去歸一
        x([1:n1],:)=x([1:n1],:)/p3;
        x([n1+1:n2],:)=x([n1+1:n2],:)/d3;
        %  反p;
        x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:);
        x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);
        x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1)*x([n1+3:n2],:);
        
        %  反d;        
        x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);
        x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1)*x([2:n1],:);
        
        %  反p;
        x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:); 
        x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x([n1+2:n2],:); 
       
        %  反d;
        x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);  
        x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:); 
        
        %  偶數(shù) 
        x([1:2:n2-1],:)=x1([1:n1],:);
        
        clear x1;
        
        %  列反變換
        
        %  歸一
        x(:,[1:n1])=x(:,[1:n1])/d3;
        x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3;
        %  反d;
        x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1);
        x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1);
        x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2)*x(:,[2:n1-1]);
        
        %  反p;        
        x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1);
        x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2)*x(:,[n1+2:n2]);
        
        %  反d;
        x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);
        x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+p1(2)*x(:,[2:n1]);
        
        %  反p;
        x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);  
        x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)*x(:,[n1+2:n2]); 
        %  奇偶
        x(:,[2:2:n2])=x(:,[n1+1:n2]);
        x(:,[1:2:n2-1])=x1(:,[1:n1]);
        clear x1;
        
        n2=n2*2; %  原大小
        n1=n2/2; %  一半大小
    end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  7.結(jié)果輸出
y=x;  %  傳輸最后結(jié)果
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  8.內(nèi)存清理
clear x;
clear flag_max;
clear layer;
clear flag_trans;
clear N;
clear n1;
clear n2;
clear s;
clear ss;
clear u;
clear d1;
clear d2;
clear d3;
clear p1;
clear p2;
clear p3;
clear sa;
clear sb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


消失矩作用的程序 
 
clear;clc;
f=50;
T=0.001;
n=1:50;
y=sin(2*pi*f*n*T);
noise=[zeros(1,17),0.2*randn(1,15),zeros(1,18)];
y_noise=y+noise;
figure(1);
subplot(2,1,1);
plot(y);
title('signal');
subplot(2,1,2);
plot(y_noise);
title('noise & signal');
[yl2,yh2]=dwt(y,'db2');
[yl10,yh10]=dwt(y,'db10');
figure(2);
subplot(2,1,1);
plot(yl2);
title('db2 low frequency signal');
subplot(2,1,2);
plot(yh2);
title('db2 high frequency signal');
figure(3);
subplot(2,1,1);
plot(yl10);
title('db10 low frequency signal');
subplot(2,1,2);
plot(yh10);
title('db10 high frequency signal');
[yl2,yh2]=dwt(y_noise,'db2');
[yl10,yh10]=dwt(y_noise,'db10');
figure(4);
subplot(2,1,1);
plot(yl2);
title('db2 low frequency noise & signal');
subplot(2,1,2);
plot(yh2);
title('db2 high frequency noise & signal');
figure(5);
subplot(2,1,1);
plot(yl10);
title('db10 low frequency noise & signal');
subplot(2,1,2);
plot(yh10);
title('db10 high frequency noise & signal'); 


小波插值與小波構造（3個程序） 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  小波構造
function casade
clear;clc;
t=3;
phi=[0,1,0];
h=wfilters('db7','r'); 
h=h*sqrt(2);
h_e=h(1,[2:2:14]);
h_o=h(1,[1:2:13]);
for m=1:15;
    stem(phi);
    drawnow;
    pause(1);
    ee=conv(h_e,phi);
    oo=conv(h_o,phi);
    phi(1,[2:2:2*length(ee)])=ee;
    phi(1,[1:2:2*length(oo)-1])=oo;
 
    
end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%  cubic_average（立方b樣條）
%  均值插值
%  初始化
s=[0 0 1 0 0]
%  正弦波
% f=50;
% ts=1/200;
% n=0:16;
% s=sin(2*pi*f*n*ts);
%  系數(shù)
se=[1/8,6/8,1/8];
so=[4/8,4/8]
%  循環(huán)
for p=1:10;
    t=length(s)-1;
    o(1:t)=s(1:t)*so(1)+s(2:t+1)*so(2);
    e(1)=s(t+1)*se(1)+s(1)*se(2)+s(2)*se(3);
    e(2:t)=s(1:t-1)*se(1)+s(2:t)*se(2)+s(3:t+1)*se(3);
    e(t+1)=s(t)*se(1)+s(t+1)*se(2)+s(1)*se(3);
    s([1:2:2*t+1])=e([1:t+1]);
    s([2:2:2*t])=o([1:t]);
    plot(s);
    drawnow;
end;
%  抽取
t=length(s); %  總長度
p=128;       %  需要點數(shù)
%  間隔
d=(t-1)/p;
%  最終尺度函數(shù)
r=s(2:d:t-1);  
%  畫圖
figure(2);
plot(r);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%  cubic_subdivision（立方插值）
%  細分插值
% %% 初始化（尺度函數(shù)）
% s=[0,0,1,0,0];
%  正弦函數(shù)
n=1:20;
f=50;
ts=1/200;
s=sin(2*pi*f*n*ts);
%  指數(shù)函數(shù)
% n=0:16;
% s=exp(n);
% %  系數(shù)
a=[-1/16,9/16,9/16,-1/16];
%  循環(huán)
for p=1:4;
    t=length(s)-1;
    o(1)=s(4)*a(1)+s(1)*a(2)+s(2)*a(3)+s(3)*a(4);
    o(2:t-1)=s(1:t-2)*a(1)+s(2:t-1)*a(2)+s(3:t)*a(3)+s(4:t+1)*a(4);
    o(t)=s(t-2)*a(4)+s(t+1)*a(3)+s(t)*a(2)+s(t-1)*a(1);
    s([1:2:2*t+1])=s([1:t+1]);
    s([2:2:2*t])=o([1:t]);
    plot(s);
    drawnow;
end;
% %  抽取
% t=length(s); %  總長度
% p=128;       %  需要點數(shù)
% 
% %  間隔
% d=(t-1)/p;
% 
% %  最終尺度函數(shù)
% r=s(2:d:t-1); 
% 
% %  畫圖
% figure(2);
% plot(r); 


小波濾波器構造和消噪程序（2個） 
 
1．重構

%  mallet_wavelet.m

%  此函數(shù)用于研究Mallet算法及濾波器設計

%  此函數(shù)僅用于消噪

 

a=pi/8;                  %角度賦初值

b=pi/8;

 

%低通重構FIR濾波器h0(n)沖激響應賦值

h0=cos(a)*cos(b);

h1=sin(a)*cos(b);

h2=-sin(a)*sin(b);

h3=cos(a)*sin(b);

 

low_construct=[h0,h1,h2,h3];

L_fre=4;                             %濾波器長度