【无标题】stft-fsst

function [sst,f,t] = fsst(x,varargin)
%FSST Fourier synchrosqueezed transform
%   SST = FSST(X) returns the Fourier synchrosqueezed transform of X. Each
%   column of SST contains the synchrosqueezed spectrum of a windowed
%   segment of X. The number of columns of SST is equal to the number of
%   samples of the input X, which is padded with zeros on each side. Each
%   spectrum is one-sided for real X and two-sided and centered for complex
%   X. By default, a Kaiser window of length 256 and a beta of 10 is used.
%   If X has fewer than 256 samples, then the Kaiser window has the same
%   length as X.
% 
%   [SST,W,S] = FSST(X) returns a vector of normalized frequencies, W,
%   corresponding to the rows of SST, and a vector of sample numbers, S,
%   corresponding to the columns of SST. Each element of S is the sample
%   number of the midpoint of a windowed segment of X.
%
%   [SST,F,T] = FSST(X,Fs) specifies the sampling frequency, Fs, in hertz,
%   as a positive scalar. The output contains sample times, T, expressed in
%   seconds and frequencies, F, expressed in hertz.
%
%   [SST,F,T] = FSST(X,Ts) specifies the sampling interval, Ts, as a 
%   <a href="matlab:help duration">duration</a>. Ts is the time between samples of X. T has the same units 
%   as the duration. The units of F are in cycles/unit time of the
%   duration.
%
%   [...] = FSST(X,...,WINDOW) when WINDOW is a vector, divides X into
%   segments of the same length as WINDOW and then multiplies each segment
%   by WINDOW.  If WINDOW is an integer, a Kaiser window with length WINDOW
%   and a beta of 10 is used. If WINDOW is not specified, the default
%   Kaiser window is used.
%
%   FSST(...) with no output arguments plots the synchrosqueezed transform
%   of the input vector x on the current figure.
%
%   FSST(...,FREQLOCATION) controls where MATLAB displays the frequency
%   axis on the plot. This string can be either 'xaxis' or 'yaxis'.
%   Setting FREQLOCATION to 'yaxis' displays frequency on the y-axis and
%   time on the x-axis.  The default is 'xaxis', which displays frequency
%   on the x-axis. FREQLOCATION is ignored when output arguments are
%   specified.
%
%   % Example 1: 
%   %   Compute and plot the Fourier synchrosqueezed transform (SST) of a
%   % signal that consists of two chirps. 
%   fs = 3000;
%   t = 0:1/fs:1-1/fs;   
%   x1 = 2*chirp(t,500,t(end),1000);
%   x2 = chirp(t,400,t(end),800); 
%   fsst(x1+x2,fs,'yaxis')
%   title('Magnitude of Fourier Synchrosqueezed Transform of Two Chirps')
%
%   % Compute and plot the short-time Fourier transform.
%   [stft,f,t] = spectrogram(x1+x2,kaiser(256,10),255,256,fs);
%   figure
%   h = imagesc(t,f,abs(stft));
%   xlabel('Time (s)') 
%   ylabel('Frequency (Hz)')
%   title('Magnitude of Short-time Fourier Transform of Two Chirps')
%   h.Parent.YDir = 'normal';
%
%   % Example 2
%   %   Compute the Fourier synchrosqueezed transform of a speech signal.
%   load mtlb
%   fsst(mtlb,Fs,hann(256),'yaxis')
%
%   See also ifsst, tfridge, spectrogram, duration.

% Copyright 2015-2016 MathWorks, Inc.

narginchk(1,4);
nargoutchk(0,3);

% Parse inputs. Fs is populated as used throughout the bulk of the code
% even if Ts is specified. fNorm specifies normalized frequencies.
[Fs,Ts,win,fNorm,freqloc] = parseInputs(x,varargin{:});

validateInputs(x,Fs,Ts,win);

% Parameters based on window size - noverlap is fixed so the transform is
% invertible.
nwin = length(win);
nfft = nwin;
noverlap = nwin-1;

% Convert to column vectors
x = signal.internal.toColIfVect(x);
win = win(:);

% Compute the output time vector (one time per sample point of the input)
if fNorm
  tout = (1:length(x));
else
  tout = (0:length(x)-1)/Fs;
end

% cast to enforce precision rules (we already checked that the inputs are
% numeric.
% cast to enforce precision rules
if (isa(x,'single') || isa(win,'single'))
  x = single(x);
  win = single(win);
  Fs = single(Fs);
  tout = single(tout);
else
  Fs = double(Fs);
end
  
% Pad the signal vector x
if mod(nwin,2)
  xp = [zeros((nwin-1)/2,1) ; x ; zeros((nwin-1)/2,1)];
else
  xp = [zeros((nwin)/2,1) ; x ; zeros((nwin-2)/2,1)];
end

nxp = length(xp);

% Place xp into columns for the STFT
xin = signal.internal.reassignSpec.getSTFTColumns(xp,nxp,nwin,noverlap,Fs);

% Compute the STFT
[sstout,fout] = computeDFT(bsxfun(@times,win,xin),nfft,Fs); 
stftc = computeDFT(bsxfun(@times,signal.internal.reassignSpec.dtwin(win,Fs),xin),nfft,Fs);

clear xin;

% Compute the reassignment vector
fcorr = -imag(stftc./ sstout);
fcorr(~isfinite(fcorr)) = 0;
fcorr = bsxfun(@plus,fout,fcorr);
tcorr = bsxfun(@plus,tout,zeros(size(fcorr)));

clear stftc;

% Mulitply STFT by a linear phase shift to produce the modified STFT
m = floor(nwin/2);
inds = 0:nfft-1;
ez = exp(-1i*2*pi*m*inds/nfft)';
sstout = bsxfun(@times,sstout,ez); 

% Reassign the modified STFT
options.range = 'twosided';
options.nfft = nfft;
sstout = signal.internal.reassignSpec.reassignSpectrum(sstout, fout, tout, fcorr, tcorr, options);

% Reduce to one-sided spectra if the input is real, otherwise return a
% two-sided (centered) spectra.
if isreal(x)
  fout = psdfreqvec('npts',nfft,'Fs',Fs,'Range','half');
  sstout = sstout(1:length(fout),:);
else
  % Centered spectra
  sstout = signal.internal.reassignSpec.centerest(sstout);
  fout = signal.internal.reassignSpec.centerfreq(fout);
end

% Scale fout and tout if the input is a duration object
if ~isempty(Ts)
  [~,~,timeScale] = getDurationandUnits(Ts);
  tout = tout*timeScale;
  fout = fout/timeScale;
end

if nargout == 0
  plotsst(tout,fout,sstout,Ts,fNorm,freqloc);
else
  sst = sstout;
  f = fout;
  t = tout(:)';
end
%--------------------------------------------------------------------------
function [Fs,Ts,win,fNorm,freqloc] = parseInputs(x,varargin)
% Set defaults
Fs = [];
Ts = [];
win = kaiser(min(256,length(x)),10); %#ok<*NASGU>
fNorm = false;
freqloc = '';

% Parse optional inputs
if nargin > 1 
  if isduration(varargin{1})
    if isempty(varargin{1})
      % Throw error is empty duration object is supplied
      error(message('signal:fsst:EmptyDuration'));
    else
      Ts = varargin{1};
      Fs = 1/seconds(Ts);
    end
  elseif ischar(varargin{1})
    freqloc = validatestring(varargin{1},{'xaxis','yaxis'});  
  else
    if ~isempty(varargin{1})
      Fs = varargin{1};
    end
  end
end

if isempty(Fs) && isempty(Ts)
  fNorm = true;
  Fs = 2*pi;
end

if nargin > 2 
  if ischar(varargin{2})
    freqloc = validatestring(varargin{2},{'xaxis','yaxis'});  
  elseif ~isempty(varargin{2})
    win = varargin{2};
    if isscalar(win)
      validateattributes(win,{'numeric'},{'positive'},'fsst','WINDOW');
      win = kaiser(double(win),10);
    end
  end
end

if nargin > 3
  freqloc = validatestring(varargin{3},{'xaxis','yaxis'}); 
end
%--------------------------------------------------------------------------
function validateInputs(x,Fs,Ts,win)

validateattributes(x,{'single','double'},...
  {'nonsparse','finite','nonnan','vector'},'fsst','X');
validateattributes(Fs,{'numeric'},...
    {'real','positive','finite','nonnan','scalar'},'fsst','Fs');
validateattributes(win,{'single','double'},...
  {'real','finite','nonnegative','nonnan','vector'},'fsst','WINDOW');  
  
if ~isempty(Ts)
  dt = getDurationandUnits(Ts);
  validateattributes(dt,{'numeric'},...
    {'real','positive','finite','nonnan','scalar'},'fsst','Ts');
end

% Check X has at least 2 samples
if length(x) < 2
  error(message('signal:fsst:MustBeMinLengthX'));
end

% Check WINDOW has at least 2 samples
if length(win) < 2
  error(message('signal:fsst:MustBeMinLengthWin'));
end

% Check window length is not more than the length of the input signal.
if length(win) > length(x)
  error(message('signal:fsst:WinLength'));
end
%--------------------------------------------------------------------------
function [dt,units,timeScale] = getDurationandUnits(Ts)
% This function returns the sampling interval and a format string
% The Units string is only for plotting.
tsformat = Ts.Format;
% Use first character of format string to determine correct
% duration object method.
tsformat = tsformat(1);
% Using the same time units as engunits. Units in engunits are
% not localized.
% time_units = {'secs','mins','hrs','days','years'};
switch tsformat
    case 's'
        dt = seconds(Ts);
        units = 'sec';
        timeScale = 1;
    case 'm'
        dt = minutes(Ts);
        units = 'min';
        timeScale = 1/seconds(minutes(1));
    case 'h'
        dt = hours(Ts);
        units = 'hr';
        timeScale = 1/seconds(hours(1));
    case 'd'
        dt = days(Ts);
        units = 'day';
        timeScale = 1/seconds(days(1));
    case 'y'
        dt = years(Ts);
        units = 'year';
        timeScale = 1/seconds(years(1));
end
%--------------------------------------------------------------------------
function plotsst(t,f,sst,Ts,fNorm,freqloc)
  % Plot the FSST in the current figure
  if fNorm
    xlbl = getString(message('signal:fsst:SampNum'));
    ylbl = [getString(message('signal:sigtools:getfrequnitstrs:NormalizedFrequency')) '  (\times\pi rad/sample)'];
    f = f/pi;
  elseif isempty(Ts) 
    % No duration object specified
    % Scale range and add corresponding label (e.g. MHz)
    [f,~,uf] = engunits(f,'unicode');
    [t,~,ut] = engunits(t,'unicode','time');
    xlbl = [getString(message('signal:spectrogram:Time')) ' (' ut ')'];
    ylbl = getfreqlbl([uf 'Hz']);
  else
    [~,timeUnit] = getDurationandUnits(Ts);
    xlbl = [getString(message('signal:spectrogram:Time')) ' (',timeUnit,'s)'];
    ylbl = [getString(message('signal:sigtools:getfrequnitstrs:Frequency')) ' (cycles/',timeUnit,')'];
  end
  
  
  if strcmp(freqloc,'yaxis')
      h = imagesc(t,f,abs(sst));
      set(h.Parent, 'Ydir', 'normal')
      xlabel(xlbl); 
      ylabel(ylbl);
  else
      h = imagesc(f,t,abs(sst)');
      set(h.Parent, 'Ydir', 'normal')
      xlabel(ylbl); 
      ylabel(xlbl);
  end
  
  title(getString(message('signal:fsst:titleFSST')));

function xrec = ifsst(sst,varargin)
%IFSST Inverse Fourier synchrosqueezed transform
%   XREC = IFSST(SST) returns the inverse of the Fourier synchrosqueezed
%   transform matrix, SST. XREC is a 1-by-N vector, where N is the number
%   of columns in SST. XREC is reconstructed using the default window and
%   the entire time-frequency plane.
%
%   XREC = IFSST(SST,WINDOW) reconstructs the signal whose SST was computed
%   by FSST using WINDOW. WINDOW is a positive vector containing the window
%   function or a scalar specifying the length of a Kaiser window with a
%   beta parameter of 10. You can leave WINDOW unspecified or as an empty
%   array in IFSST and the default Kaiser window will be used. 
%
%   XREC = IFSST(SST,WINDOW,F,FREQRANGE) uses the frequency vector F and
%   the two-element vector FREQRANGE to invert the synchrosqueezed
%   transform. F is the vector of frequencies corresponding to the rows of
%   SST. The synchrosqueezed transform is inverted for the bins in SST with
%   frequencies in FREQRANGE. FREQRANGE is a two-element vector with
%   strictly increasing values and must lie in the range of F.
%
%   XREC = IFSST(SST,WINDOW,IRIDGE) inverts the synchrosqueezed transform
%   along the time-frequency ridges specified by the index vector or matrix
%   IRIDGE. IRIDGE is an output of TFRIDGE. If IRIDGE is a matrix, IFSST
%   first inverts the synchrosqueezed transform along the first column of
%   IRIDGE, then iteratively reconstructs along subsequent columns of
%   IRIDGE. XREC has the same size as IRIDGE.
%
%   XREC = IFSST(SST,WINDOW,IRIDGE,'NumFrequencyBins',NBINS) specifies the
%   number of frequency bins on each side of IRIDGE to use in the
%   reconstruction. If the index of the time-frequency ridge +/- NBINS
%   exceeds the range of frequency bins at any time step, IFSST truncates
%   the reconstruction at the first or last frequency bin. If unspecified,
%   NBINS defaults to 4. The name-value pair 'NumFrequencyBins', NBINS is
%   valid only if you provide the IRIDGE input.
%
%   % Example 1
%   %   Extract and reconstruct each of two signal components using the
%   %   Fourier synchrosqueezed transform.  
%   fs = 3000; 
%   t = 0:1/fs:1-1/fs; 
%   x1 = 2*chirp(t,500,t(end),1000); 
%   x2 = chirp(t,400,t(end),800); 
%   [sst,f] = fsst(x1+x2,fs);
%
%   %   Extract a ridge for each component.
%   [~,iridge] = tfridge(sst,f,'NumRidges',2);
%
%   %   Reconstruct a waveform for each ridge.
%   xrec = ifsst(sst,[],iridge);
%
%   %   Plot the extracted waveforms.
%   figure, subplot(2,1,1)
%   plot(t,xrec(:,1))
%   title('Waveform Extraction')
%   xlabel('Time (s)'), ylabel('Amplitude')
%   legend('Chirp 1')
%   grid
%   subplot(2,1,2)
%   plot(t,xrec(:,2),'Color',[0.85 0.325 0.098])
%   xlabel('Time (s)'), ylabel('Amplitude')
%   legend('Chirp 2')
%   grid
%
%   % Example 2
%   %   Compute the Fourier synchrosqueezed transform of a speech signal.
%   %   Invert the transform to reconstruct the signal, plot the
%   %   reconstructed signal, and compare the original and reconstructed 
%   %   signals using the L-infinity norm.
%   
%   %   Load the speech signal and compute the synchrosqueezed transform.
%   load mtlb
%   [sst,f] = fsst(mtlb,Fs);
%
%   %   Invert the synchrosqueezed transform and reconstruct the signal.
%   xrec = ifsst(sst);
%
%   %   Plot the orignal and reconstructed signals.
%   t = (0:length(mtlb)-1)/Fs;
%   plot(t,mtlb,t,xrec)
%   xlabel('Time (s)')
%   legend('Original','Reconstructed')
%   title('Original and Reconstructed Speech Signals')
%   
%   %   Compute the L-infinity norm of the difference between original and
%   %   reconstructed signals.
%   Linf = norm(abs(mtlb-xrec),Inf)
%
%   See also fsst, tfridge.

% Copyright 2015-2016 MathWorks, Inc.

narginchk(1,5);
nargoutchk(0,1);

[win,iridge,f,freqrange,nfft,nbins,fReal] = parseInputs(sst,varargin{:});

validateInputs(sst,win,iridge,f,freqrange,nbins);

% Cast to enforce precision rules
% The output xrec will be single only if the input sst is single. If the
% window used was single, sst would be single.
win = double(win);

% Convert to column vectors
win = win(:);
f = f(:);
freqrange = freqrange(:);

% Compute window scale factor, frequency interval, and real input flag
g0 = win(round(length(win)/2)); 
df = 1/nfft;

% Determine the last index to use in the second side of the spectrum.
% For an even length fft, the one-sided spectrum has (nfft/2)+1 elements,
% while for an odd length fft, the one-sided spectrum has (nfft+1)/2
% elements.
if mod(nfft,2)
  iend = size(sst,1);
else
  iend = size(sst,1)-1;
end
% Extract ridges if provided
if ~isempty(iridge)
  xrec = zeros(size(iridge));
  for j = 1:size(iridge,2)
    Mask = zeros(size(sst));
    for i = 1:size(iridge,1)
       Mask(max(iridge(i,j)-nbins,1):min(iridge(i,j)+nbins,size(sst,1)),i) = 1;
    end
    % Do the inverse 
    maskSST = sst.*Mask;
    if fReal
      xrec(:,j) = 1/g0*df*real(sum(maskSST)+ ...
        conj(sum(maskSST(2:iend,:))))';
    else
      xrec(:,j) = 1/g0*df*(sum(maskSST)).';  
    end
  end
elseif ~isempty(freqrange);
  [~,i1] = min(abs(f - freqrange(1)));  
  [~,i2] = min(abs(f - freqrange(2)));
   % Do the inverse 
    if fReal
      xrec = 1/g0*df*real(sum(sst(i1:i2,:)) + ...
        conj(sum(sst(max(i1,2):min(i2,iend),:))));
    else
      xrec = 1/g0*df*(sum(sst(i1:i2,:)));
    end
  xrec = xrec(:);
else
  % Do the inverse 
  if fReal
    xrec = 1/g0*df*real(sum(sst) + conj(sum(sst(2:iend,:))));
  else
    xrec = 1/g0*df*(sum(sst));
  end
  xrec = xrec(:);
end

%--------------------------------------------------------------------------
function [win,iridge,f,freqrange,nfft,nbins,fReal] = parseInputs(sst,varargin)

sstSize=size(sst);

if sstSize(2) < 256
   win = kaiser(sstSize(2),10);
else
   win = kaiser(256,10);
end

iridge = [];
f = [];
freqrange = [];
nbins = 4; %#ok<*NASGU>

if nargin > 1
  % The first optional argument is WINDOW
    win = getWin(varargin{1},win);
end

if nargin == 3
  % This is a ridge
    iridge = varargin{2};
end

if nargin == 4
  % This F and FREQRANGE
  if ~ischar(varargin{3})
    f = varargin{2};
    freqrange = varargin{3};
  else
    error(message('signal:ifsst:NVMustBeEven'));
  end
end

if nargin == 5
  iridge = varargin{2};
  % This is ridge and name-value pair
  validatestring(varargin{3},{'NumFrequencyBins'});
  nbins = varargin{4};
end

nfft = length(win);

if sstSize(1) == nfft
  fReal = false;
else
  fReal = true;
end

%--------------------------------------------------------------------------
function win = getWin(newwin,currentWin)
  if isempty(newwin)
    win = currentWin;
  elseif isscalar(newwin)
      validateattributes(newwin,{'single','double'},{'positive'},'ifsst','WINDOW');
      win = kaiser(newwin,10);
  else
    win = newwin;
  end
%--------------------------------------------------------------------------
function validateInputs(sst,win,iridge,f,freqrange,nbins)

validateattributes(sst,{'single','double'},...
  {'nonsparse','finite','nonnan','nonempty'},'ifsst','SST');
validateattributes(win,{'single','double'},...
  {'real','finite','nonnegative','nonnan','vector'},'ifsst','WINDOW');

if ~isempty(iridge)
  validateattributes(iridge,{'single','double'},...
    {'real','positive','integer','nonsparse','finite','nonnan'},'ifsst','IRIDGE');
end
if ~isempty(f) || ~isempty(freqrange)
  validateattributes(f,{'single','double'},...
    {'real','finite','nonnan','vector'},'ifsst','F');
  validateattributes(freqrange,{'single','double'},...
    {'real','finite','nonnan','vector'},'ifsst','FREQRANGE');
end
validateattributes(nbins,{'numeric'},...
  {'real','positive','integer','finite','nonnan','scalar','nonempty'},'ifsst','BINS');

% check that sst is a matrix. fsst has a minimum signal size of 2, which
% will produce a 2x2 array.
if isvector(sst)
  error(message('signal:ifsst:SSTisVector'));
end  

% Window length matches the size of sst
sizeSST = size(sst);
winLen = length(win);

% Expected window length for a real signal
if mod(winLen,2)
  winLenReal = 2*sizeSST(1)-1; 
else
  winLenReal = 2*(sizeSST(1)-1); 
end

if ~(winLen == sizeSST(1) || winLen == winLenReal)
  error(message('signal:ifsst:winLen'));
end

% Size of f matches SST
if ~isempty(f) && ~(length(f) == size(sst,1))
  error(message('signal:ifsst:fSize'));
  % The length of F must match the number of rows of SST. 
end

% Size of iridge matches SST
if ~isempty(iridge) && ~(size(iridge,1) == size(sst,2))
  error(message('signal:ifsst:iridgeSize'));
  % The number of rows of IRIDGE must match the number of rows of SST.
end

% Length of freqrange is 2, values are within the range of f and strictly
% increasing
if ~isempty(freqrange) && (~(length(freqrange) == 2) || ... 
  freqrange(2) <= freqrange(1) || freqrange(1) < f(1) || freqrange(2) > f(end)) 
  error(message('signal:ifsst:freqrange'));
end