www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/eml/imodwpt.m

    function xrec = imodwpt(cfs,varargin)
%MATLAB Code Generation Private Function

%   Copyright 1995-2016 The MathWorks, Inc.
%#codegen

narginchk(1,4);
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
% Input must a real-valued matrix with no Infs or NaNs
validateattributes(cfs,{'double'},{'real','finite'},'imodwpt','CFS');
% The terminal level must be at least j=1, which means that cfs must have
% two rows
coder.internal.assert(~isvector(cfs),'Wavelet:modwt:InvalidCFSSize');
NumPackets = coder.internal.indexInt(size(cfs,1));
% The number of rows in the matrix must be a power of two
coder.internal.assert(coder.internal.sizeIsPow2(NumPackets), ...
    'Wavelet:modwt:InvalidTermSize');
%The level of the transform
level = coder.internal.indexInt(log2(double(NumPackets)));
%Parse inputs
defaultWav = 'fk18';
[~,Lo,Hi] = parseModwptInputs(level,defaultWav,varargin{:});
Lo = Lo/sqrt(2);
Hi = Hi/sqrt(2);
% If the coefficient length is less than the filter length, need to
% periodize the signal in order to use the DFT algorithm
N = coder.internal.indexInt(size(cfs,2));
%For the edge case where the number of samples is less than the scaling
%filter
nLo = coder.internal.indexInt(numel(Lo));
if N < nLo
    ncopies = nLo - N;
    Nrep = (1 + ncopies)*N;
else
    ncopies = coder.internal.indexInt(0);
    Nrep = N;
end
% Allocate storage for a possibly expanded data array.
cfs1 = coder.nullcopy(zeros(size(cfs,1),Nrep));
% Copy one page of elements.
cfs1(:,1:N) = cfs;
if ncopies > 0
    for k = 1:ncopies
        offset = k*N;
        for j = 1:N
            cfs1(:,offset + j) = cfs1(:,j);
        end
    end
end
% Obtain the DFT of the filters
G = fft(Lo,double(Nrep));
H = fft(Hi,double(Nrep));
mcfs = coder.internal.indexInt(size(cfs1,1));
upcfs = zeros(mcfs/2,Nrep);
for jj = level:-1:1
    kk = ONE;
    index = ZERO;
    twopjjm1 = eml_lshift(ONE,jj - 1); 
    for nn = 0:twopjjm1-1
        index = index + 1;
        if eml_bitand(nn,ONE) == 0
            upcfs(index,:) = EvenInvert(cfs1(kk,:),cfs1(kk+1,:),G,H,jj);
        else
            upcfs(index,:) = OddInvert(cfs1(kk,:),cfs1(kk+1,:),G,H,jj);
        end
        kk = kk + 2;
    end
    mcfs = mcfs/2;
    cfs1(1:mcfs,:) = upcfs(1:mcfs,:);
end
%Ensure output length matches the number of columns in the input
xrec = cfs1(1,1:N);

%--------------------------------------------------------------------------

function evencfs = EvenInvert(even,odd,G,H,J)
N = coder.internal.indexInt(numel(even));
ONE = coder.internal.indexInt(1);
evendft = fft(even);
odddft = fft(odd);
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
for k = 1:N
    idx = 1 + mod(upfactor*(k - 1),N);
    Gupk = conj(G(idx));
    Hupk = conj(H(idx));
    % Reuse storage for evendft: evendft = Gup.*evendft + Hup.*odddft.
    evendft(k) = Gupk*evendft(k) + Hupk*odddft(k);
end
evencfs = real(ifft(evendft));

%--------------------------------------------------------------------------

function oddcfs = OddInvert(even,odd,G,H,J)
N = coder.internal.indexInt(numel(even));
ONE = coder.internal.indexInt(1);
evendft = fft(even);
odddft = fft(odd);
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
for k = 1:N
    idx = 1 + mod(upfactor*(k - 1),N);
    Gupk = conj(G(idx));
    Hupk = conj(H(idx));
    % Reuse storage for odddft: odddft = Gup.*odddft + Hup.*evendft.
    odddft(k) = Gupk*odddft(k) + Hupk*evendft(k);
end
oddcfs = real(ifft(odddft));

%--------------------------------------------------------------------------