Matlab 实现小波变换

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('haar')%Haar 变换的分解和重构滤波器的长度为2waveinfo('haar');[phi,psi,xval]=wavefun('haar',10);xaxis=zeros(size(xval));subplot(121);plot(xval,phi,'k',xval,xaxis,'--k');axis([0 1 -1.5 1.5]);axis square;title('haar scaling function');subplot(122);plot(xval,psi,'k',xval,xaxis,'--k');axis([0 1 -1.5 1.5]);axis square;title('haar wavelet function');

f=magic(4)%使用Haar 的单尺度小波变换[c1,s1]=wavedec2(f,1,'haar')[c2,s2]=wavedec2(f,2,'haar')

function [varargout]=wavefilter(wname,type)
%Check the input and output argumentserror (nargchk(1,2,nargin));if (nargin==1&nargout~=4)|(nargin==2&nargout~=2)    error('Invalid number of output arguments');endif nargin == 1 & ~ischar(wname)    error('WNAME must be a string');end
if nargin == 2 & ~ischar(type)    error('TYPE must be a string');end
switch lower(wname)    case {'haar','db1'}        ld=[1 1]/sqrt(2); hd=[-1 1]/sqrt(2);        lr=ld; hr=-hd;    case 'db4'        ld=[-1.059740178499728e-002 3.288301166698295e-002 ...            3.084138183598697e-002 -1.870348117188811e-001 ...            -2.798376941698385e-002 6.308807679295904e-001 ...            7.14846575525415e-001 2.303778133088552e-001];        t=(0:7);        hd=ld; hd(end:-1:1)=cos(pi*t).*ld;        lr=ld; lr(end:-1:1)=ld;        hr=cos(pi*t).*ld;    case 'sym4'        ld=[-7.576571478927333e-002 -2.963552764599851e-002 ...            4.976186676320155e-001 8.03738758059161e-001 ...            2.978577956052774e-001 -9.921954357684722e-002 ...            -1.260396726203783e-002 3.222310060404270e-002];        t=(0:7);        hd=ld; hd(end:-1:1)=cos(pi*t).*ld;        lr=ld; lr(end:-1:1)=ld;        hr=cos(pi*t).*ld;    case 'bior6.8'        ld=[0 1.908831736481291e-003 -1.914286129088767e-003 ...            -1.699063986760234e-002 1.193456527972926e-002 ...            4.973490349094079e-002 -7.726317316720414e-002 ...            -9.405920349573646e-002 4.207962846098268e-001 ...            8.259229974584023e-001 4.207962846098268e-001 ...            -9.405920349573646e-002 -7.726317316720414e-002 ...            4.973490349094079e-002 1.193456527972926e-002 ...            -1.699063986760234e-002 -1.914286129088767e-003 ...            1.908831736481291e-003];        hd=[0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...            -7.872200106262882e-002 4.036797903033992e-002 ...            4.178491091502746e-001 -7.589077294536542e-002 ...            4.178491091502746e-001 4.036797903033992e-002 ...            -7.872200106262882e-002 -1.446750489679015e-002 ...            1.442628250562444e-002 0 0 0 0];        t=(0:17);        lr=cos(pi*(t+1)).*hd;        hr=cos(pi*t).*ld;    case 'jpep9.7'        ld=[0 0.02674875741080976 -0.01686411844287495 ...            -0.07822326652898785 0.266861184428732 ...            0.6029490182363579 0.266861184428732 ...            -0.07822326652898785 -0.01686411844287495 ...            0.02674875741080976];        hd=[0 -0.09127176311424948 0.05754352622849957 ...            0.5912717631142470 -1.115087052456994 ...            0.5912717631142470 0.05754352622849957 ...            -0.09127176311424948 0 0];        t=(0:9);        lr=cos(pi*(t+1)).*hd;        hr=cos(pi*t).*ld;    otherwise        error('Unrecognizable wavelet name (WNAME).');end
%Output the requested filtersif(nargin == 1)    varargout(1:4)={ld,hd,lr,hr};else    switch lower(type(1))        case 'd'        varargout={ld,hd};        case 'r'        varargout={lr,hr};        otherwise            error('Unrecognizable filter TYPE.')     endend

function [c,s]=wavefast(x,n,varargin)
error(nargchk(3,4,nargin));if nargin == 3    if ischar(varargin{1})        [lp,hp]=wavefilter(varargin{1},'d');    else        error('Missing wavelet name');    endelse    lp=varargin{1};hp=varargin{2};endf1=length(lp);sx=size(x);if (ndims(x)~=2)|(min(sx<2))|~isreal(x)|~isnumeric(x)    error('X must be a real, numeric matrix');end
if (ndims(lp)~=2)|~isreal(lp)|~isnumeric(lp) ...    |(ndims(hp)~=2)|~isreal(hp)|~isnumeric(hp) ...    |(f1~=length(hp)) | rem(f1,2)~=0    error(['LP and HP must be even and equal length real ,'...'numeric filter vectors']);end
if ~isreal(n)|~isnumeric(n)|(n<1)|(n>log2(max(sx)))    error(['N must be a real scalar between 1 and' ...'log2(max(size(x)))']);endc=[];s=sx;app=double(x);
%for each decomposition...for i=1:n    [app,keep]=symextend(app,f1);
    rows=symconv(app,hp,'row',f1,keep);    coefs=symconv(rows,hp,'col',f1,keep);    c=[coefs(:)' c];    s=[size(coefs);s];    coefs=symconv(rows,lp,'col',f1,keep);    c=[coefs(:)' c];    rows=symconv(app,lp,'row',f1,keep);    coefs=symconv(rows,hp,'col',f1,keep);    c=[coefs(:)' c];    app=symconv(rows,lp,'col',f1,keep);endc=[app(:)' c];s=[size(app);s];
function [y,keep]=symextend(x,f1)keep=floor((f1+size(x)-1)/2);y=padarray(x,[(f1-1) (f1-1)],'symmetric','both');

function [ratio, maxdiff]=fwtcompare(f,n,wname)tic;[c1,s1]=wavedec2(f,n,wname);reftime=toc;tic;[c2,s2]=wavefast(f,n,wname);t2=toc;

%compare the resultsratio=t2/(reftime+eps);maxdiff=abs(max(c1-c2)); 

function y=symconv(x,h,type,f1,keep)if strcmp(type,'row')    y=conv2(x,h);    y=y(:,1:2:end);    y=y(:,f1/2+1:f1/2+keep(2));else    y=conv2(x,h');    y=y(1:2:end,:);    y=y(f1/2+1:f1/2+keep(1),:);end 
%比较两种方法下的速度clear allclcf=imread('D:\图像库\DIP3E_CH07_Original_Images\Fig0701.tif','tif');[ratio,maxdifference]=fwtcompare(f,5,'db4')

clear allclcf=imread('D:\图像库\DIP3E_CH07_Original_Images\Fig0701.tif');imshow(f);[c,s]=wavefast(f,1,'sym4');figure;wave2gray(c,s,-6);[nc,y]=wavecut('a',c,s);figure;wave2gray(nc,s,-6);edges=abs(waveback(nc,s,'sym4'));figure;imshow(mat2gray(edges));

%子程序function w=wave2gray(c,s,scale,border)error(nargchk(2,4,nargin));if(ndims(c)~=2)|(size(c,1)~=1)    error('c must be a row vector');endif(ndims(s)~=2)|~isreal(s)|~isnumeric(s)|(size(s,2)~=2)    error('s must be a real,numeric two-column array');endelements=prod(s,2);
if (length(c)<elements(end))|...        ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))    error(['[c s] must be a standard wavelet' ...'decomposition structure']);endif(nargin>2)&(~isreal(scale)|~isnumeric(scale))    error('scale must be a real,numeric scalar');end
if(nargin>3)&(~ischar(border))    error('border must be character string');end
if nargin==2    scale=1;%default scaleendif nargin<4    border='absord';%default borderend
absflag=scale<0;scale=abs(scale);if scale==0    scale=1;end
[cd,w]=wavecut('a',c,s);w=mat2gray(w);cdx=max(abs(cd(:)))/scale;if absflag    cd=mat2gray(abs(cd),[0,cdx]);    fill=0;else    cd=mat2gray(cd,[-cdx,cdx]);    fill=0.5;end%build gray image one decomposition at a timefor i=size(s,1)-2:-1:1    ws=size(w);    h=wavecopy('h',cd,s,i);    pad=ws-size(h);    frontporch=round(pad/2);    h=padarray(h,frontporch,fill,'pre');    h=padarray(h,pad-frontporch,fill,'post');    v=wavecopy('v',cd,s,i);    pad=ws-size(v);    frontporch=round(pad/2);    v=padarray(v,frontporch,fill,'pre');    h=padarray(v,pad-frontporch,fill,'post');    d=wavecopy('d',cd,s,i);    pad=ws-size(d);    frontporch=round(pad/2);    d=padarray(d,frontporch,fill,'pre');    d=padarray(d,pad-frontporch,fill,'post');
    %Add 1 pixel white border    switch lower(border)        case 'append'            w=padarray(w,[1 1],1,'post');            h=padarray(h,[1 0],1,'post');            v=padarray(v,[0 1],1,'post');        case 'absord'            w(:,end)=1;            w(end,:)=1;            h(end,:)=1;            v(:,end)=1;        otherwise            error('Unrecognized border parameter');    end    w=[w h;v d];end
if nargout==0    imshow(w);%display resultend

function [varargout]=waveback(c,s,varargin)
error(nargchk(3,5,nargin));error(nargchk(1,2,nargout));
if (ndims(c)~=2)|(size(c,1)~=1)    error('c must be a row vector');endif(ndims(s)~=2)|~isreal(s)|~isnumeric(s)|(size(s,2)~=2)    error('s must be a real,numeric two-column array');endelements=prod(s,2);
if (length(c)<elements(end))|...~(elements(1)+3*sum(elements(2:end-1))>=elements(end))     error(['[c s] must be a standard wavelet' ...    'decomposition structure']);end
nmax=size(s,1)-2;wname=varargin{1};filterchk=0;nchk=0;
switch nargin    case 3        if ischar(wname)            [lp,hp]=wavefilter(wname,'r');            n=nmax;        else            error('updefined filter');        end        if nargout~=1            error('wrong number of output arguments');        end    case 4        if ischar(wname)            [lp,hp]=wavefilter(wname,'r');            n=varargin{2};            nchk=1;        else            lp=varargin{1};            hp=varargin{2};            filterchk=1;            n=nmax;            if nargout~=1                error('wrong number of output arguments');            end        endcase 5    lp=varargin{1};    hp=varargin{2};    filterchk=1;    n=varargin{3};    nchk=1;    otherwise        error('improper number fo input arguments');end
f1=length(lp);if filterchk    if (ndims(lp)~=2)|~isreal(lp)|~isnumeric(lp) ...        |(ndims(hp)~=2)|~isreal(hp)|~isnumeric(hp) ...        |(f1~=length(hp)) | rem(f1,2)~=0    error(['LP and HP must be even and equal length real ,'...        'numeric filter vectors']);    endend
if nchk & (~isnumeric(n)|~isreal(n))    error('n must be a real numeric');endif (n>nmax)|(n<1)    error('Invalid number (N) of reconstructions requested');endif (n~=nmax) & (nargout~=2)    error('Not enough output arguments');end
nc=c;ns=s;nnmax=nmax;for i=1:n    a=symconvup(wavecopy('a',nc,ns),lp,lp,f1,ns(3,:))+ ...        symconvup(wavecopy('h',nc,ns,nnmax),hp,lp,f1,ns(3,:))+ ...        symconvup(wavecopy('v',nc,ns,nnmax),lp,hp,f1,ns(3,:))+ ...        symconvup(wavecopy('d',nc,ns,nnmax),hp,hp,f1,ns(3,:));
        nc=nc(4*prod(ns(1,:))+1:end);        nc=[a(:)' nc];        ns=ns(3:end,:);        ns=[ns(1,:);ns];        nnmax=size(ns,1)-2;end
if nargout==1;    a=nc;    nc=repmat(0,ns(1,:));    nc(:)=a;endvarargout{1}=nc;if nargout==2    varargout{2}=ns;end
function z=symconvup(x,f1,f2,fln,keep)
y=zeros([2 1].*size(x));y(1:2:end,:)=x;y=conv2(y,f1');z=zeros([1 2].*size(y));z(:,1:2:end)=y;z=conv2(z,f2);z=z(fln-1:fln+keep(1)-2,fln-1:fln+keep(2)-2);
function [nc,y]=wavecut(type,c,s,n)error(nargchk(3,4,nargin));if nargin==4    [nc,y]=wavework('cut',type,c,s,n);else    [nc,y]=wavework('cut',type,c,s);End
function y=wavecopy(type,c,s,n)error(nargchk(3,4,nargin));if nargin==4    [nc,y]=wavework('cut',type,c,s,n);else    [nc,y]=wavework('cut',type,c,s);End
function y=wavepaste(type,c,s,n)error(nargchk(5,5,nargin))nc=wavework('paste',type,c,s,n,x);
function [varargout]=wavework(opcode,type,c,s,n,x)error(nargchk(4,6,nargin));if (ndims(c)~=2)|(size(c,1)~=1)    error('c must be a row vector');endif(ndims(s)~=2)|~isreal(s)|~isnumeric(s)|(size(s,2)~=2)    error('s must be a real,numeric two-column array');end
elements=prod(s,2);
if (length(c)<elements(end))|...        ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))    error(['[c s] must be a standard wavelet' ...        'decomposition structure']);end
if strcmp (lower(opcode(1:3)),'pas') & nargin < 6    error('Not enough input arguments');end
if nargin<5    n=1;end

nmax=size(s,1)-2;

aflag=(lower(type(1))=='a');if ~aflag & (n>nmax)    error('N exceeds the decompositon in [c,s]');end
switch lower(type(1))    case 'a'        nindex=1;        start=1;        stop=elements(1);        ntst=nmax;    case {'h','v','d'}        switch type            case 'h'                offset=0;            case 'v'                offset=1;            case 'd'                offset=2;        end        nindex=size(s,1)-n;        start=elements(1)+3*sum(elements(2:nmax-n+1))+ ...            offset*elements(nindex)-1;        stop=start+elements(nindex)-1;        ntst=n;    otherwise        error('Type must begin with "a","h","v",or "d"');end
switch lower(opcode)    case {'copy','cut'}        y=repmat(0,s(nindex,:));        y(:)=c(start:stop);        nc=c;        if strcmp(lower(opcode(1:3)),'cut')            nc(start:stop)=0;            varargout={nc,y};        else            varargout={y};        end    case 'paste'        if prod(size(x))~=elements(end-ntst)            error('x is not sized for the requested paste');        else            nc=c;            nc(start:stop)=x(:);            varargout={nc};        end    otherwise        error('unrecognized OPCODE');end