www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/eml/modwtmra.m
function mra = modwtmra(w,varargin) %MATLAB Code Generation Library Function % Copyright 1995-2016 The MathWorks, Inc. %#codegen narginchk(1,4); coder.internal.assert(~isvector(w),'Wavelet:modwt:MRASize'); validateattributes(w,{'double'},{'real','nonnan','finite'},mfilename); ncw = coder.internal.indexInt(size(w,2)); % get the size of the output coefficients cfslength = coder.internal.indexInt(ncw); J0 = coder.internal.indexInt(size(w,1)) - 1; nullinput = zeros(1,cfslength); N = cfslength; ZERO = coder.internal.indexInt(0); ONE = coder.internal.indexInt(1); %Parse input arguments defaultLev = ZERO; defaultWave = 'sym4'; [~,Lo,Hi,reflection] = parseModwtInputs( ... defaultLev,defaultWave,varargin{:}); % Scale the scaling and wavelet filters for the MODWT Lo = Lo./sqrt(2); Hi = Hi./sqrt(2); % Adjust final output length if MODWT obtained with 'reflection' if reflection N = eml_rshift(N,ONE); end nLo = coder.internal.indexInt(numel(Lo)); if cfslength < nLo ncopies = nLo - cfslength; cfslength = (1 + ncopies)*cfslength; else ncopies = ZERO; end % Allocate storage for a possibly expanded data array. ww = coder.nullcopy(zeros(size(w,1),cfslength)); % Copy one page of elements. ww(:,1:ncw) = w; if ncopies > 0 for k = 1:ncopies offset = k*ncw; for j = 1:ncw ww(:,offset + j) = ww(:,j); end end end % Obtain the DFT of the filters G = fft(Lo,double(cfslength)); H = fft(Hi,double(cfslength)); % Allocate array for MRA mra = zeros(J0+1,N); for J = J0:-1:1 details = imodwtDetails(w(J,:),nullinput,J,G,H,cfslength); mra(J,:) = details(1:N); end % scalingcoefs = w(J0+1,:); smooth = imodwtSmooth(w(J0+1,:),nullinput,G,H,cfslength,J0); mra(J0+1,:) = smooth(1:N); %----------------------------------------------------------------------- function details = imodwtDetails(coefs,nullinput,lev,Lo,Hi,N) v = nullinput; w = coefs; for jj = lev:-1:1 vout = imodwtrec(v,w,Lo,Hi,jj); v = vout; w = nullinput; end details = v(1:N); %----------------------------------------------------------------------- function smooth = imodwtSmooth(scalingcoefs,nullinput,Lo,Hi,N,J0) v = scalingcoefs; for J = J0:-1:1 vout = imodwtrec(v,nullinput,Lo,Hi,J); v = vout; end smooth = v(1:N); %---------------------------------------------------------------------- function Vout = imodwtrec(Vin,Win,G,H,J) N = coder.internal.indexInt(length(Vin)); ONE = coder.internal.indexInt(1); upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1); Vhat = fft(Vin); What = fft(Win); for k = 1:N idx = 1 + mod(upfactor*(k - 1),N); Gupk = conj(G(idx)); Hupk = conj(H(idx)); % Reuse storage for Vhat: Vhat = Gup.*Vhat + Hup.*What Vhat(k) = Gupk*Vhat(k) + Hupk*What(k); end Vout = real(ifft(Vhat)); %--------------------------------------------------------------------------