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