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

    function [W,packetlevels,F] = modwptdetails(x,varargin)
%MATLAB Code Generation Library Function

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

narginchk(1,6);
ONE = coder.internal.indexInt(1);
N = coder.internal.indexInt(numel(x));
defaultLev = coder.internal.indexInt(min(4,floor(log2(double(N)))));
defaultWav = 'fk18';
[J,Lo,Hi,fulltree] = parseModwptInputs( ...
    defaultLev,defaultWav,varargin{:});
validateinputs(x,J,Lo,Hi,fulltree);
% This is the initial DFT length for the input
Lo = Lo/sqrt(2);
Hi = Hi/sqrt(2);
% If the signal length is less than the filter length, need to
% periodize the signal in order to use the DFT algorithm
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.
xx = coder.nullcopy(zeros(1,Nrep));
xx(1:N) = x;
if ncopies > 0
    % x = [x,repmat(x,1,ncopies)];
    for k = 1:ncopies
        offset = k*N;
        for j = 1:N
            xx(offset + j) = xx(j);
        end
    end
end
% Obtain the DFT of the filters
% G is the DFT of the scaling filter, H is DFT of wavelet filter
G = fft(Lo,double(Nrep));
H = fft(Hi,double(Nrep));
% Obtain DFT of original data
X = fft(xx,double(Nrep));
% Create array to hold wavelet packets and packet levels
% Initially create full tree
mcfs = eml_lshift(ONE,J+1) - 2;
cfs = zeros(mcfs,Nrep,'like',X);
cfs(1,:) = X;
% MODWPT algorithm
p2 = ONE;
for kk = 1:J
    nIdx = p2;
    p2 = p2*2; % p2 = 2^kk
    % Determine first packet for a given level
    jj = p2 - 1;
    if kk > 1
        Idx2 = lengthPacketLevels(kk - 1);
        Idx = Idx2 + 1 - nIdx;
    else
        Idx = ONE;
    end
    nnhi = p2/2 - 1;
    for nn = 0:nnhi
        X = cfs(Idx + nn,:);
        [vhat,what] = modwptdecxcorr(X,G,H,kk);
        if eml_bitand(nn,ONE) == ONE % odd
            cfs(jj + 2*nn,:) = what;
            cfs(jj + 2*nn + 1,:) = vhat;
        else % even
            cfs(jj + 2*nn + 1,:) = what;
            cfs(jj + 2*nn,:) = vhat;
        end
    end
end
% Now p2 = 2^J, but let's be explicit.
pow2J = eml_lshift(ONE,J);
W1 = coder.nullcopy(zeros(mcfs,N));
for k = 1:mcfs
    v = real(ifft(cfs(k,:)));
    % Ensure output length matches signal length
    W1(k,:) = v(1:N);
end
if fulltree
    W = W1;
else
    offset = mcfs - pow2J;
    W = coder.nullcopy(zeros(pow2J,N,'like',W1));
    for j = 1:N
        for i = 1:pow2J
            W(i,j) = W1(offset + i,j);
        end
    end
end
if nargout > 1
    [packetlevels,F] = calcModwptPLandF(J,fulltree);
end

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

function [Vhat,What] = modwptdecxcorr(X,G,H,J)
% [Vhat,What] = modwtfft(X,G,H,J)
ONE = coder.internal.indexInt(1);
N = coder.internal.indexInt(length(X));
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
Vhat = coder.nullcopy(X);
What = coder.nullcopy(X);
for k = 1:N
    idx = 1 + mod(upfactor*(k - 1),N);
    Vhat(k) = abs(G(idx)).^2*X(k);
    What(k) = abs(H(idx)).^2*X(k);
end

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

function validateinputs(x,J,Lo,Hi,fulltree)
coder.inline('always');
coder.internal.prefer_const(J,Lo,Hi,fulltree);
%Input must be real-valued, double with no Infs or NaNs
validateattributes(x,{'double'},{'real','finite'},'modwpt','X');
%Input must be 1-D
coder.internal.assert(isvector(x), ...
    'Wavelet:modwt:OneD_Input');
%Input must contain at least two samples
N = numel(x);
coder.internal.assert(N >= 2,'Wavelet:modwt:LenTwo');
% J is the transform level
% validateattributes(J,{'double'},{'integer','positive'},'modwpt','LEVEL');
%Check the transform level does not exceed the maximum
coder.internal.assert(J <= floor(log2(N)),'Wavelet:modwt:MRALevel');

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