极大重叠离散小波变换 分解 matlab代码

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