wavebac.m

%%wavebac源程序

function [varargout]=waveback(c,s,varargin)
%waveback函数执行一个二维n阶水平部分或完整的小波分解结构[c,s]的重构

%syntax:
%y=waveback(c,s,'wname');输出FWT的逆矩阵y
%y=waveback(c,s,lr,hr);使用低通和高通重构滤波器或名为'wname'的波形滤波器
%[nc,ns]=waveback(c,s,'wname',n);输出新的小波
%[nc,ns]=waveback(c,s,lr,hr,n);n步重构后的分解结构[nc,ns]

%检查输入和输出参数
error(nargchk(3,5,narargin));
error(nargchk(1,2,nargout));

if(ndims(c)~=2) | (size(c,1)~=1)
    error('C must be a row vector.');
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

%在[c,s]中的最大水平
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('Undefined 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
        end 
    case 5
        lp=varargin{1};hp=varargin{2};filterchk=1;
        n=varargin{3};nchk=1;
    otherwise
        error('Improper number of 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.']);
    end 
end

if nchk & (~isnumeric(n) | isreal(n))
    error('N must be a real numeric.');
end
if (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,:);
        nnmax=size(ns,1)-2;
end

%为完成重构，重定输出格式为二维
if nargout==1
    a=nc;nc=repmat(0,ns(1,:));nc(:)=a;
end 

varargout{1}=nc;
if nargout==2
    varargout{2}=ns;
end