function w = md_dec(x,varargin)
%MODWT Maximal overlap discrete wavelet transform.
narginchk(1,5);
% Validate that data is real, 1-D double-precision
% with no NaNs or Infs
validateattributes(x,{'double'},{'real','nonnan','finite'});
%Input must be 1-D
if (~isrow(x) && ~iscolumn(x))
error(message('Wavelet:modwt:OneD_Input'));
end
%Input must contain at least two samples
if (numel(x)<2)
error(message('Wavelet:modwt:LenTwo'));
end
% Convert data to row vector
x = x(:)';
% Record original data length
datalength = length(x);
%Parse input arguments
params = parseinputs(datalength,varargin{:});
%Check that the level of the transform does not exceed floor(log2(numel(x))
J = params.J;
Jmax = floor(log2(datalength));
if (J <= 0) || (J > Jmax) || (J ~= fix(J))
error(message('Wavelet:modwt:MRALevel'));
end
boundary = params.boundary;
if (~isempty(boundary) && ~strcmpi(boundary,'reflection'))
error(message('Wavelet:modwt:Invalid_Boundary'));
end
% increase signal length if 'reflection' is specified
if strcmpi(boundary,'reflection')
x = [x flip(x)];
end
% obtain new signal length if needed
siglen = length(x);
Nrep = siglen;
% If wavelet specified as a string, ensure that wavelet is orthogonal
if (isfield(params,'wname') && ~isfield(params,'Lo'))
[~,~,Lo,Hi] = wfilters(params.wname);
wtype = wavemngr('type',params.wname);
if (wtype ~= 1)
error(message('Wavelet:modwt:Orth_Filt'));
end
end
%If scaling and wavelet filters are specified as vectors, ensure they
%satisfy the orthogonality conditions
if (isfield(params,'Lo') && ~isfield(params,'wname'))
filtercheck = wavelet.internal.CheckFilter(params.Lo,params.Hi);
if ~filtercheck
error(message('Wavelet:modwt:Orth_Filt'));
end
Lo = params.Lo;
Hi = params.Hi;
end
% Scale the scaling and wavelet filters for the MODWT
Lo = Lo./sqrt(2);
Hi = Hi./sqrt(2);
% Ensure Lo and Hi are row vectors
Lo = Lo(:)';
Hi = Hi(:)';
% If the signal length is less than the filter length, need to
% periodize the signal in order to use the DFT algorithm
if (siglen<numel(Lo))
x = [x repmat(x,1,ceil(numel(Lo)-siglen))];
Nrep = numel(x);
end
% Allocate coefficient array
w = zeros(J+1,Nrep);
% Obtain the DFT of the filters
G = fft(Lo,Nrep);
H = fft(Hi,Nrep);
%Obtain the DFT of the data
Vhat = fft(x);
% Main MODWT algorithm
for jj = 1:J
[Vhat,What] = modwtdec(Vhat,G,H,jj);
w(jj,:) = ifft(What);
end
w(J+1,:) = ifft(Vhat);
% Truncate data to length of boundary condition
w = w(:,1:siglen);
%----------------------------------------------------------------------
function [Vhat,What] = modwtdec(X,G,H,J)
% [Vhat,What] = modwtfft(X,G,H,J)
N = length(X);
upfactor = 2^(J-1);
Gup = G(1+mod(upfactor*(0:N-1),N));
Hup = H(1+mod(upfactor*(0:N-1),N));
Vhat = Gup.*X;
What = Hup.*X;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function params = parseinputs(siglen,varargin)
% % Parse varargin and check for valid inputs
% % First convert any strings to char arrays
% [varargin{:}] = convertStringsToChars(varargin{:});
% % Assign defaults
% params.boundary = [];
% params.J = floor(log2(siglen));
% params.wname = 'sym4';
%
% % Check for 'reflection' boundary
% tfbound = strcmpi(varargin,'reflection');
%
% % Determine if 'reflection' boundary is specified
% if any(tfbound)
% params.boundary = varargin{tfbound>0};
% varargin(tfbound>0) = [];
% end
%
% % If boundary is the only input in addition to the data, return with
% % defaults
% if isempty(varargin)
% return;
% end
%
% % Only remaining char variable must be wavelet name
% tfchar = cellfun(@ischar,varargin);
% if (nnz(tfchar) == 1)
% params.wname = varargin{tfchar>0};
% end
%
% % Only scalar input must be the level
% tfscalar = cellfun(@isscalar,varargin);
%
% % Check for numeric inputs
% tffilters = cellfun(@isnumeric,varargin);
%
% % At most 3 numeric inputs are supported
% if nnz(tffilters)>3
% error(message('Wavelet:modwt:Invalid_Numeric'));
% end
%
% % There's one numeric argument and it's not a wavelet level
% if (nnz(tffilters)==1) && (nnz(tfscalar) == 0)
% error(message('Wavelet:FunctionInput:InvalidLoHiFilters'));
% end
%
% % If there are at least two numeric inputs, the first two must be the
% % scaling and wavelet filters
% if (nnz(tffilters)>1)
% idxFilt = find(tffilters,2,'first');
% params.Lo = varargin{idxFilt(1)};
% params.Hi = varargin{idxFilt(2)};
% params = rmfield(params,'wname');
%
% if (length(params.Lo) < 2 || length(params.Hi) < 2)
% error(message('Wavelet:modwt:Invalid_Filt_Length'));
% end
%
% end
%
% % Any scalar input must be the level
% if any(tfscalar)
% params.J = varargin{tfscalar>0};
% end
%
% % If the user specifies a filter, use that instead of default wavelet
% if (isfield(params,'Lo') && any(tfchar))
% error(message('Wavelet:FunctionInput:InvalidWavFilter'));
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params = parseinputs(siglen,varargin)
% Parse varargin and check for valid inputs
% First convert any strings to char arrays
% Assign defaults
params.boundary = [];
params.J = floor(log2(siglen));
params.wname = 'sym4';
% Check for 'reflection' boundary
tfbound = strcmpi(varargin,'reflection');
% Determine if 'reflection' boundary is specified
if any(tfbound)
params.boundary = varargin{tfbound>0};
varargin(tfbound>0) = [];
end
% If boundary is the only input in addition to the data, return with
% defaults
if isempty(varargin)
return;
end
% Only remaining char variable must be wavelet name
tfchar = cellfun(@ischar,varargin);
if (nnz(tfchar) == 1)
params.wname = varargin{tfchar>0};
end
% Only scalar input must be the level
tfscalar = cellfun(@isscalar,varargin);
% Check for numeric inputs
tffilters = cellfun(@isnumeric,varargin);
% At most 3 numeric inputs are supported
if nnz(tffilters)>3
error(message('Wavelet:modwt:Invalid_Numeric'));
end
% There's one numeric argument and it's not a wavelet level
if (nnz(tffilters)==1) && (nnz(tfscalar) == 0)
error(message('Wavelet:FunctionInput:InvalidLoHiFilters'));
end
% If there are at least two numeric inputs, the first two must be the
% scaling and wavelet filters
if (nnz(tffilters)>1)
idxFilt = find(tffilters,2,'first');
params.Lo = varargin{idxFilt(1)};
params.Hi = varargin{idxFilt(2)};
params = rmfield(params,'wname');
if (length(params.Lo) < 2 || length(params.Hi) < 2)
error(message('Wavelet:modwt:Invalid_Filt_Length'));
end
end
% Any scalar input must be the level
if any(tfscalar)
params.J = varargin{tfscalar>0};
end
% If the user specifies a filter, use that instead of default wavelet
if (isfield(params,'Lo') && any(tfchar))
error(message('Wavelet:FunctionInput:InvalidWavFilter'));
end
- 1
- 2
- 3
- 4
- 5
- 6
前往页