www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/wsst.m

    function [sst,f]  = wsst(x,varargin)
%Wavelet Synchrosqueezed Transform
%   SST = wsst(X) returns the wavelet synchrosqueezed transform for the 1-D
%   real-valued signal, X. X must have at least 4 samples. The
%   synchrosqueezed transform uses 32 voices per octave and the number of
%   octaves is floor(log2(numel(X)))-1. The transform uses the analytic
%   Morlet wavelet by default. SST is a Na-by-N matrix where Na is the
%   number of scales, 32*(floor(log2(numel(X)))-1), and N is the number of
%   samples in X.
%
%   [SST,F] = wsst(X) returns a vector of frequencies, F, in cycles/sample
%   corresponding to the rows of SST.
%
%   [...] = wsst(X,Fs) specifies the sampling frequency, Fs, in hertz as a
%   positive scalar. If you specify the sampling frequency, WSST returns
%   the frequencies in hertz.
%
%   [...] = wsst(X,Ts) uses the positive <a href="matlab:help duration">duration</a>, Ts, to compute the
%   scale-to-frequency conversion, F. Ts is the time between samples of X.
%   F has units of cycles/unit time where the unit of time is the same time
%   unit as the duration.
%
%   [...] = wsst(...,WAV) uses the analytic wavelet specified by WAV to
%   compute the synchrosqueezed transform. Valid choices for WAV are
%   'amor' and 'bump' for the analytic Morlet and bump wavelet. If
%   unspecified, WAV defaults to 'amor'.
%
%   [...] = wsst(...,'VoicesPerOctave',NV) specifies the number of voices
%   per octave as a positive even integer between 10 and 48. The number of
%   scales is the product of the number of voices per octave and the number
%   of octaves. If unspecified, NV defaults to 32 voices per octave. You
%   can specify the 'VoicesPerOctave' name-value pair anywhere in the input
%   argument list after the signal X.
%
%   [...] = wsst(...,'ExtendSignal',EXTENDFLAG) specifies whether to
%   symmetrically extend the signal by reflection to mitigate boundary
%   effects. EXTENDFLAG can be one of the following options [ true |
%   {false}]. If unspecified, EXTENDSIGNAL defaults to false.  You can
%   specify the 'ExtendSignal' name-value pair anywhere in the input
%   argument list after the signal X.
%
%   wsst(...) with no output arguments plots the wavelet synchrosqueezed
%   transform as a function of time and frequency. If you do not specify a
%   sampling frequency or interval, the synchrosqueezed transform is
%   plotted in cycles/sample. If you supply a sampling frequency, Fs, the
%   synchrosqueezed transform is plotted in hertz. If you supply a <a href="matlab:help duration">duration</a>
%   as a sampling interval, the synchrosqueezed transform is plotted
%   in cycles/unit time where the time unit is the same as the duration.
%
%   % Example 1:
%   %   Obtain the wavelet synchrosqueezed transform of a quadratic chirp.
%   %   The chirp is sampled at 1000 Hz.
%   load quadchirp;
%   [sst,f] = wsst(quadchirp,1000);
%   hp = pcolor(tquad,f,abs(sst));
%   hp.EdgeColor = 'none';
%   title('Wavelet Synchrosqueezed Transform');
%   xlabel('Time'); ylabel('Hz');
%
%   % Example 2:
%   %   Obtain the wavelet synchrosqueezed transform of the sunspot
%   %   data. Specify the sampling interval to be 1 for one sample per
%   %   year.
%   load sunspot.dat;
%   wsst(sunspot(:,2),years(1))
%
%   See also iwsst, wsstridge, duration
narginchk(1,8);
nbSamp = numel(x);
x = x(:)';
validateattributes(x,{'double'},{'row','finite','real'},'wsst','X');
if numel(x)<4
    error(message('Wavelet:synchrosqueezed:NumInputSamples'));
end
params = parseinputs(nbSamp,varargin{:});
nv = params.nv;
noct = params.noct;
% Create scale vector
na = noct*params.nv;


% If sampling frequency is specified, dt = 1/fs
if (isempty(params.fs) && isempty(params.Ts))
    % The default is 1 for normalized frequency
    dt = params.dt;
    Units = '';
elseif (~isempty(params.fs) && isempty(params.Ts))
    % Accept the sampling frequency in hertz
    fs = params.fs;
    dt = 1/fs;
    Units = '';
elseif (isempty(params.fs) && ~isempty(params.Ts))
    % Get the dt and Units from the duration object
    [dt,Units] = getDurationandUnits(params.Ts);
    
    
end

a0 = 2^(1/nv);
scales = a0.^(1:na);
NbSc = numel(scales);

% Construct time series to analyze, pad if necessary
meanSIG = mean(x);
x = x - meanSIG;
NumExten = 0;

if params.pad
    %Pad the time series symmetrically
    np2 = nextpow2(nbSamp);
    NumExten = 2^np2-nbSamp;
    x = wextend('1d','symw',x,NumExten,'b');
end

%Record data length plus any extension
N = numel(x);

%Create frequency vector for CWT computation
omega = (1:fix(N/2));
omega = omega.*((2.*pi)/N);
omega = [0., omega, -omega(fix((N-1)/2):-1:1)];

% Compute FFT of the (padded) time series
xdft = fft(x);
[psift,dpsift]  = sstwaveft(params.WAV,omega,scales,params.wavparam);

%Obtain CWT coefficients and derivative
cwtcfs = ifft(repmat(xdft,NbSc,1).*psift,[],2);
dcwtcfs = ifft(repmat(xdft,NbSc,1).*dpsift,[],2);

%Remove padding if any
cwtcfs = cwtcfs(:,NumExten+1:end-NumExten);
dcwtcfs = dcwtcfs(:,NumExten+1:end-NumExten);

%Compute the phase transform
phasetf = imag(dcwtcfs./cwtcfs)./(2*pi);


% Threshold for synchrosqueezing
phasetf(abs(phasetf)<params.thr) = NaN;

% Create frequency vector for output
log2Nyquist = log2(1/(2*dt));
log2Fund = log2(1/(nbSamp*dt));
freq = 2.^linspace(log2Fund,log2Nyquist,na);

Tx = 1/nv*sstalgo(cwtcfs,phasetf,params.thr);


if (nargout == 0)
    
    plotsst(Tx,freq,dt,params.engunitflag,params.normalizedfreq,Units);
else
    sst = Tx;
    f = freq;
end

%-------------------------------------------------------------------------
function [wft,dwft] = sstwaveft(WAV,omega,scales,wavparam)
%   Admissible wavelets are:
%    - MORLET wavelet (A) - 'morl':
%        PSI_HAT(s) = exp(-(s-s0).^2/2) * (s>0)
%        Parameter: s0, default s0 = 6.
%   - Bump wavelet:  'bump':
%       PSI_HAT(s) = exp(1-(1/((s-mu)^2./sigma^2))).*(abs((s-mu)/sigma)<1)
%       Parameters: mu,sigma.
%       default:    mu=5, sigma = 0.6.
%   Normalized to have unit magnitude at the peak frequency of the wavelet

NbSc = numel(scales);
NbFrq = numel(omega);
wft = zeros(NbSc,NbFrq);

switch WAV
    case 'amor'
        
        cf = wavparam;
        
        for jj = 1:NbSc
            expnt = -(scales(jj).*omega - cf).^2/2.*(omega > 0);
            wft(jj,:) = exp(expnt).*(omega > 0);
        end
        
    case 'bump'
        
        mu = wavparam(1);
        sigma = wavparam(2);
        
        
        for jj = 1:NbSc
            w = (scales(jj)*omega-mu)./sigma;
            expnt = -1./(1-w.^2);
            daughter = exp(1)*exp(expnt).*(abs(w)<1-eps(1));
            daughter(isnan(daughter)) = 0;
            wft(jj,:) = daughter;
        end
        
end

%Compute derivative
omegaMatrix = repmat(omega,NbSc,1);
dwft = 1j*omegaMatrix.*wft;

%-------------------------------------------------------------------------
function plotsst(Tx,F,dt,engunitflag,isfreqnormalized,Units)

if ~isempty(Units)
    freqUnits = Units(1:end-1);    
end

t = 0:dt:(size(Tx,2)*dt)-dt;
if engunitflag && isfreqnormalized
    frequnitstrs = wgetfrequnitstrs;
    freqlbl = frequnitstrs{1};
    xlbl = 'Samples';
elseif engunitflag && ~isfreqnormalized
    [F,~,uf] = wengunits(F,'unicode');
    freqlbl = wgetfreqlbl([uf 'Hz']);
    [t,~,ut] = wengunits(t,'unicode','time');
    xlbl = [getString(message('Wavelet:getfrequnitstrs:Time')) ' (' ut ')'];
    
else
    freqlbl = getString(message('Wavelet:synchrosqueezed:FreqLabel'));
    freqlbl = ...
        [freqlbl '/' freqUnits ')'];
    xlbl = getString(message('Wavelet:synchrosqueezed:Time'));
    xlbl = [xlbl ' (' Units ')'];
end


h = pcolor(t,F,abs(Tx));
h.EdgeColor = 'none';
shading interp;
ylabel(freqlbl); xlabel(xlbl);
title(getString(message('Wavelet:synchrosqueezed:SynchrosqueezedTitle')));



%-------------------------------------------------------------------------
function params = parseinputs(nbSamp,varargin)
% Set defaults.
params.fs = [];
params.dt = 1;
params.Ts = [];
params.sampinterval = false;
params.engunitflag = true;
params.WAV = 'amor';
params.wavparam = 6;
params.thr = 1e-8;
params.nv = 32;
params.noct = floor(log2(nbSamp))-1;
params.pad = false;
params.normalizedfreq = true;

% Error out if there are any calendar duration objects
tfcalendarDuration = cellfun(@iscalendarduration,varargin);
if any(tfcalendarDuration)
    error(message('Wavelet:FunctionInput:CalendarDurationSupport'));
end

tfsampinterval = cellfun(@isduration,varargin);

if (any(tfsampinterval) && nnz(tfsampinterval) == 1)
    params.sampinterval = true;
    params.Ts = varargin{tfsampinterval>0};
    if (numel(params.Ts) ~= 1 ) || params.Ts <= 0 || isempty(params.Ts)
        error(message('Wavelet:FunctionInput:PositiveScalarDuration'));
    end
    
    params.engunitflag = false;
    params.normalizedfreq = false;
    varargin(tfsampinterval) = [];
end

%Look for Name-Value pairs
numvoices = find(strncmpi('voicesperoctave',varargin,1));

if any(numvoices)
    params.nv = varargin{numvoices+1};
    %validate the value is logical
    validateattributes(params.nv,{'numeric'},{'positive','scalar',...
        'even','>=',10,'<=',48},'wsst','VoicesPerOctave');
    varargin(numvoices:numvoices+1) = [];
    if isempty(varargin)
        return;
    end
end


extendsignal = find(strncmpi('extendsignal',varargin,1));

if any(extendsignal)
    params.pad = varargin{extendsignal+1};
    
    if ~isequal(params.pad,logical(params.pad))
        error(message('Wavelet:FunctionInput:Logical'));
    end
    varargin(extendsignal:extendsignal+1) = [];
    if isempty(varargin)
        return;
    end
end


% Only scalar left must be sampling frequency or sampling interval
% Only scalar left must be sampling frequency
tfsampfreq = cellfun(@(x) (isscalar(x) && isnumeric(x)),varargin);

if (any(tfsampfreq) && (nnz(tfsampfreq) == 1) && ~params.sampinterval)
    params.fs = varargin{tfsampfreq};
    validateattributes(params.fs,{'numeric'},{'positive'},'wsst','Fs');
    params.normalizedfreq = false;
    params.engunits = true;
elseif any(tfsampfreq) && params.sampinterval
    error(message('Wavelet:FunctionInput:SamplingIntervalOrDuration'));
elseif nnz(tfsampfreq)>1
    error(message('Wavelet:FunctionInput:Invalid_ScalNum'));
end

%Only char variable left must be wavelet
tfwav = cellfun(@ischar,varargin);
if (nnz(tfwav) == 1)
    params.WAV = varargin{tfwav>0};
    params.WAV = validatestring(params.WAV,{'bump','amor'},'wsst','WAV');
elseif nnz(tfwav)>1
    error(message('Wavelet:FunctionInput:InvalidChar'));
    
end

if strncmpi(params.WAV,'bump',1)
    params.wavparam = [5 1];
end



%------------------------------------------------------------------------
function Tx = sstalgo(cwtcfs,phasetf,gamma)

M = size(cwtcfs,1);
N = size(cwtcfs,2);
log2Fund = log2(1/N);
log2Nyquist = log2(1/2);
iRow = real(1 + floor(M/(log2Nyquist-log2Fund)*(log2(phasetf)-log2Fund)));
idxphasetf = find(iRow>0 & iRow<=M & ~isnan(iRow));
idxcwtcfs = find(abs(cwtcfs)>gamma);
idx = intersect(idxphasetf,idxcwtcfs);
iCol = repmat(1:N,M,1);
Tx = accumarray([iRow(idx) iCol(idx)],cwtcfs(idx),size(cwtcfs));