www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/icwt.m
function xrec = icwt(wt,varargin) %ICWT Inverse continuous 1-D wavelet transform % XREC = ICWT(WT) inverts the continuous wavelet transform (CWT) % coefficient matrix, WT. ICWT assumes that you used the default Morse % (3,60) wavelet and default scales in obtaining the CWT. % % XREC = ICWT(WT,WNAME) uses the wavelet WNAME in the inversion. WNAME is a % valid wavelet name: 'morse', 'amor', or 'bump'. The inverse CWT must use % the same wavelet used in the CWT. % % XREC = ICWT(WT,F,FREQRANGE) inverts the CWT over the two-element % frequency range specified in FREQRANGE. The elements of FREQRANGE must be % stricly increasing and contained in the range of the frequency vector F. % F is assumed to be the scale-to-frequency conversion obtained in CWT. % % XREC = ICWT(WT,PERIOD,PERIODRANGE) inverts the CWT over the two-element % range of periods in PERIODRANGE. The elements of PERIODRANGE must be % strictly increasing and contained in the range of the period vector % PERIOD. PERIOD is an array of durations obtained from CWT with a duration % input. % % XREC = ICWT(...,'TimeBandwidth',TB) use the positive scalar % time-bandwidth parameter, TB, to invert the CWT using the Morse wavelet. % The symmetry parameter (gamma) of the Morse wavelet is assumed to equal % 3. The inverse CWT must use the same time-bandwidth value used in the % CWT. % % XREC = ICWT(...,'WaveletParameters',PARAM) uses the parameters PARAM to % specify the Morse wavelet used in the inversion of the CWT. PARAM is a % two-element vector. The first element is the symmetry parameter (gamma) % and the second parameter is the time-bandwidth parameter. The inverse CWT % must use the same wavelet parameters used in the CWT. % % XREC = ICWT(...,'SignalMean',MEAN) adds the scalar MEAN to the output of % ICWT. Because the continuous wavelet transform does not preserve the % signal mean, the inverse CWT is a zero-mean signal by default. Note that % adding a non-zero MEAN to a frequency- or period-limited reconstruction % adds a zero-frequency component to the reconstruction. % % XREC = ICWT(...,'VoicesPerOctave',NV) specifies the number of voices per % octave used in obtaining the CWT. If you input a frequency vector or % array of durations, you cannot specify the VoicesPerOctave name-value % pair. The number of voices per octave is determined by the frequency or % duration vector. If you do not specify the number of voices per octave or % a frequency or duration vector, ICWT uses the default of 10. NV is an % even integer between 4 and 48 and must agree with the value used in % obtaining the CWT. % % % %Example 1: % % Obtain the CWT of the Kobe earthquake data. Invert the CWT % % and compare the result with the original signal. % load kobe; % sigmean = mean(kobe); % wt = cwt(kobe); % xrec = icwt(wt,'SignalMean',sigmean); % plot((1:numel(kobe))./60,kobe); % xlabel('mins'); ylabel('nm/s^2'); % hold on; % plot((1:numel(kobe))./60,xrec,'r'); % legend('Inverse CWT','Original Signal'); % % %Example 2: % % Reconstruct a frequency-localized approximation to the Kobe % % earthquake data by extracting information from the CWT % % corresponding to frequencies in the range of [0.030, 0.070] Hz. % load kobe; % [wt,f] = cwt(kobe,1); % xrec = icwt(wt,f,[0.030 0.070],'SignalMean',mean(kobe)); % subplot(211) % plot(kobe); grid on; % title('Original Data'); % subplot(212) % plot(xrec); grid on; % title('Bandpass Filtered Reconstruction [0.030 0.070] Hz'); % % See also CWT narginchk(1,8); nargoutchk(0,1); validateattributes(wt,{'numeric'},{'nonempty','finite'}); if isrow(wt) || iscolumn(wt) error(message('Wavelet:cwt:InvalidCWTSize')); end %Get the number of scales Na = size(wt,1); params = parseinputs(Na,varargin{:}); ds = params.ds; if ~isempty(params.f) ds = getScType(1./params.f); idxbegin = find(params.f>=params.freqrange(2),1,'last'); idxend = find(params.f<=params.freqrange(1),1,'first'); idxZero = setdiff(1:Na,idxbegin:idxend); wt(idxZero,:) = 0; elseif ~isempty(params.periods) ds = getScType(params.periods); idxbegin = find(params.periods<=params.periodrange(1),1,'last'); idxend = find(params.periods>=params.periodrange(2),1,'first'); idxZero = setdiff(1:Na,idxbegin:idxend); wt(idxZero,:) = 0; end morseparams = [params.ga params.be]; cpsi = wavelet.internal.admConstant(params.wavname,morseparams); a0 = 2^ds; Wr = 2*log(a0)*(1/cpsi)*real(wt); xrec = sum(Wr,1); xrec = xrec+params.SignalMean; %--------------------------------------------------------------------- function ds = getScType(scales) DF2 = mean(diff(log(scales)/log(scales(1)),2)); if abs(DF2) < sqrt(eps) ds = mean(diff(log2(scales))); else error(message('Wavelet:cwt:UnsupportedScales')); end %----------------------------------------------------------------------- function params = parseinputs(Na,varargin) % Set defaults. params.wavname = 'morse'; params.ga = 3; params.be = 20; params.SignalMean = 0; params.ds = 1/10; params.freqrange = []; params.f = []; params.periodrange = []; params.periods = []; params.duration = false; morseparams = find(strncmpi(varargin,'waveletparameters',1)); timeBandwidth = find(strncmpi(varargin,'timebandwidth',1)); if any(morseparams) && any(timeBandwidth) error(message('Wavelet:cwt:paramsTB')); end if (any(morseparams) && (nnz(morseparams) == 1)) morseParameter = varargin{morseparams+1}; validateattributes(morseParameter,{'numeric'},{'numel',2,... 'positive','nonempty'},'icwt','MorseParameters'); varargin(morseparams:morseparams+1) = []; params.ga = morseParameter(1); tb = morseParameter(2); validateattributes(params.ga,{'numeric'},{'scalar',... 'positive','>=',1},'icwt','gamma'); validateattributes(tb,{'numeric'},{'scalar',... '>',params.ga,'<=',40*params.ga},... 'icwt','TimeBandWidth'); params.be = tb/params.ga; end if (any(timeBandwidth) && (nnz(timeBandwidth) == 1)) params.timebandwidth = varargin{timeBandwidth+1}; validateattributes(params.timebandwidth,{'numeric'},{'scalar',... 'positive','>' 3, '<',120},'icwt','TimeBandwidth'); params.ga = 3; params.be = params.timebandwidth/params.ga; if params.be>40 error(message('Wavelet:cwt:TBupperbound')); end varargin(timeBandwidth:timeBandwidth+1) = []; end tfmean = find(strncmpi(varargin,'SignalMean',1)); if (any(tfmean) && nnz(tfmean) == 1) params.SignalMean = varargin{tfmean+1}; validateattributes(params.SignalMean,{'numeric'},{'scalar','real',... 'finite'},'icwt','SignalMean'); varargin(tfmean:tfmean+1) = []; end tfduration = cellfun(@isduration,varargin); if nnz(tfduration) == 1 error(message('Wavelet:cwt:PeriodAlone')); end if (any(tfduration) && nnz(tfduration) == 2) durationidx = find(tfduration); params.periods = varargin{durationidx(1)}; params.periodrange = varargin{durationidx(2)}; [params.periods,UnitsP] = getDurationandUnits(params.periods); [params.periodrange,UnitsPR] = getDurationandUnits(params.periodrange); if ~strcmp(UnitsPR,UnitsP) error(message('Wavelet:cwt:PeriodUnits')); end validateattributes(params.periods,{'numeric'},{'positive','nonempty'... 'size',[Na 1]},'icwt','PERIOD'); psmall = params.periods(1); plarge = params.periods(end); % Check that the periodrange is contained in periods validateattributes(params.periodrange,{'numeric'},... {'positive','nonempty','finite','increasing','numel',2,... '>=',psmall,'<=',plarge},'icwt','PERIODRANGE'); end tfnumvoices = find(strncmpi(varargin,'voicesperoctave',1)); if any(tfnumvoices) && nnz(tfnumvoices)==1 nv = varargin{tfnumvoices+1}; validateattributes(nv,{'numeric'},{'positive','scalar',... 'even','>=',10,'<=',48},'icwt','VoicesPerOctave'); varargin(tfnumvoices:tfnumvoices+1) = []; params.ds = 1/nv; if isempty(varargin) return; end end tfnumeric = cellfun(@isnumeric,varargin); if nnz(tfnumeric) == 1 error(message('Wavelet:cwt:FrequencyAlone')); end if (any(tfnumeric) && nnz(tfnumeric) == 2) && ~any(tfduration) % If there are two numeric inputs, the first must be F and % the second must be FREQRANGE -- you cannot specify both % frequencies and periods idxvector = find(tfnumeric); params.f = varargin{idxvector(1)}; params.freqrange = varargin{idxvector(2)}; validateattributes(params.f,{'numeric'},{'positive','decreasing',... 'size',[Na 1]},'icwt','F'); % Find the limits of F to test against FREQRANGE fhigh = params.f(1); flow = params.f(end); validateattributes(params.freqrange,{'numeric'},... {'numel',2,'<=',fhigh,'>=',flow},'icwt','FREQRANGE'); elseif (any(tfnumeric) && nnz(tfnumeric) == 2) && any(tfduration) error(message('Wavelet:FunctionInput:FrequencyOrPeriod')); end if any(tfnumvoices) && (~isempty(params.f) || ~isempty(params.periods)) error(message('Wavelet:cwt:NVandFrequencies')); end %Only char variable left must be wavelet tfwav = cellfun(@ischar,varargin); if (nnz(tfwav) == 1) params.wavname = varargin{tfwav>0}; params.wavname = ... validatestring(params.wavname,{'morse','bump','amor'},'icwt','WAVNAME'); if ~strcmp(params.wavname,'morse') params.ga = []; params.be = []; end elseif nnz(tfwav)>1 error(message('Wavelet:FunctionInput:InvalidChar')); end