基于matlab的广义S变换仿真

function [st,t,f] = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)

% Returns the Stockwell Transform of the timeseries.

% Code by Robert Glenn Stockwell.

% DO NOT DISTRIBUTE

% BETA TEST ONLY

% Reference is "Localization of the Complex Spectrum: The S Transform"

% from IEEE Transactions on Signal Processing, vol. 44., number 4, April 1996, pages 998-1001.

%

%-------Inputs Needed------------------------------------------------

%  

%   *****All frequencies in (cycles/(time unit))!******

% "timeseries" - vector of data to be transformed

%-------Optional Inputs ------------------------------------------------

%

%"minfreq" is the minimum frequency in the ST result(Default=0)

%"maxfreq" is the maximum frequency in the ST result (Default=Nyquist)

%"samplingrate" is the time interval between samples (Default=1)

%"freqsamplingrate" is the frequency-sampling interval you desire in the ST result (Default=1)

%Passing a negative number will give the default ex.  [s,t,f] = st(data,-1,-1,2,2)

%-------Outputs Returned------------------------------------------------

%

% st     -a complex matrix containing the Stockwell transform. 

%  The rows of STOutput are the frequencies and the 

%         columns are the time values ie each column is 

%         the "local spectrum" for that point in time

%  t      - a vector containing the sampled times

%  f      - a vector containing the sampled frequencies

%--------Additional details-----------------------

%   %  There are several parameters immediately below that

%  the user may change. They are:

%[verbose]    if true prints out informational messages throughout the function.

%[removeedge] if true, removes a least squares fit parabola

%                and puts a 5% hanning taper on the edges of the time series.

%                This is usually a good idea.

%[analytic_signal]  if the timeseries is real-valued

%                      this takes the analytic signal and STs it.

%                      This is almost always a good idea.

%[factor]     the width factor of the localizing gaussian

%                ie, a sinusoid of period 10 seconds has a 

%                gaussian window of width factor*10 seconds.

%                I usually use factor=1, but sometimes factor = 3

%                to get better frequency resolution.

%   Copyright (c) by Bob Stockwell

%   $Revision: 1.2 $  $Date: 1997/07/08  $

% This is the S transform wrapper that holds default values for the function.

TRUE = 1;

FALSE = 0;

%%% DEFAULT PARAMETERS  [change these for your particular application]

verbose = TRUE;          

removeedge= FALSE;

analytic_signal =  FALSE;

factor = 1;

%%% END of DEFAULT PARAMETERS

%%%START OF INPUT VARIABLE CHECK

% First:  make sure it is a valid time_series 

%         If not, return the help message

if verbose disp(' '),end  % i like a line left blank

if nargin == 0 

   if verbose disp('No parameters inputted.'),end

   st_help

   t=0;,st=-1;,f=0;

   return

end

% Change to column vector

if size(timeseries,2) > size(timeseries,1)

timeseries=timeseries';

end

% Make sure it is a 1-dimensional array

if size(timeseries,2) > 1

   error('Please enter a *vector* of data, not matrix')

return

elseif (size(timeseries)==[1 1]) == 1

error('Please enter a *vector* of data, not a scalar')

return

end

% use defaults for input variables

if nargin == 1

   minfreq = 0;

   maxfreq = fix(length(timeseries)/2);

   samplingrate=1;

   freqsamplingrate=1;

elseif nargin==2

   maxfreq = fix(length(timeseries)/2);

   samplingrate=1;

   freqsamplingrate=1;

   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

elseif nargin==3 

   samplingrate=1;

   freqsamplingrate=1;

   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

elseif nargin==4   

   freqsamplingrate=1;

   [ minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

elseif nargin == 5

      [ minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries);

else      

   if verbose disp('Error in input arguments: using defaults'),end

   minfreq = 0;

   maxfreq = fix(length(timeseries)/2);

   samplingrate=1;

   freqsamplingrate=1;

end

if verbose 

   disp(sprintf('Minfreq = %d',minfreq))

   disp(sprintf('Maxfreq = %d',maxfreq))

   disp(sprintf('Sampling Rate (time   domain) = %d',samplingrate))

   disp(sprintf('Sampling Rate (freq.  domain) = %d',freqsamplingrate))

   disp(sprintf('The length of the timeseries is %d points',length(timeseries)))

   disp(' ')

end

%END OF INPUT VARIABLE CHECK

% If you want to "hardwire" minfreq & maxfreq & samplingrate & freqsamplingrate do it here

% calculate the sampled time and frequency values from the two sampling rates

t = (0:length(timeseries)-1)*samplingrate;

spe_nelements =ceil((maxfreq - minfreq+1)/freqsamplingrate)   ;

f = (minfreq + [0:spe_nelements-1]*freqsamplingrate)/(samplingrate*length(timeseries));

if verbose disp(sprintf('The number of frequency voices is %d',spe_nelements)),end

% The actual S Transform function is here:

st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor); 

% this function is below, thus nicely encapsulated

%WRITE switch statement on nargout

% if 0 then plot amplitude spectrum

if nargout==0 

   if verbose disp('Plotting pseudocolor image'),end

   pcolor(t,f,abs(st))

end

return

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

function st = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor); 

% Returns the Stockwell Transform, STOutput, of the time-series

% Code by R.G. Stockwell.

% Reference is "Localization of the Complex Spectrum: The S Transform"

% from IEEE Transactions on Signal Processing, vol. 44., number 4,

% April 1996, pages 998-1001.

%

%-------Inputs Returned------------------------------------------------

%         - are all taken care of in the wrapper function above

%

%-------Outputs Returned------------------------------------------------

%

% ST    -a complex matrix containing the Stockwell transform.

%  The rows of STOutput are the frequencies and the

%  columns are the time values

%

%

%-----------------------------------------------------------------------

% Compute the length of the data.

n=length(timeseries);

original = timeseries;

if removeedge

    if verbose disp('Removing trend with polynomial fit'),end

   ind = [0:n-1]';

    r = polyfit(ind,timeseries,2);

    fit = polyval(r,ind) ;

 timeseries = timeseries - fit;

    if verbose disp('Removing edges with 5% hanning taper'),end

    sh_len = floor(length(timeseries)/10);

    wn = hanning(sh_len);

    if(sh_len==0)

       sh_len=length(timeseries);

       wn = 1&[1:sh_len];

    end

    % make sure wn is a column vector, because timeseries is

   if size(wn,2) > size(wn,1)

      wn=wn';

   end

   

   timeseries(1:floor(sh_len/2),1) = timeseries(1:floor(sh_len/2),1).*wn(1:floor(sh_len/2),1);

timeseries(length(timeseries)-floor(sh_len/2):n,1) = timeseries(length(timeseries)-floor(sh_len/2):n,1).*wn(sh_len-floor(sh_len/2):sh_len,1);

  

end

% If vector is real, do the analytic signal 

if analytic_signal

   if verbose disp('Calculating analytic signal (using Hilbert transform)'),end

   % this version of the hilbert transform is different than hilbert.m

   %  This is correct!

   ts_spe = fft(real(timeseries));

   h = [1; 2*ones(fix((n-1)/2),1); ones(1-rem(n,2),1); zeros(fix((n-1)/2),1)];

   ts_spe(:) = ts_spe.*h(:);

   timeseries = ifft(ts_spe);

end  

% Compute FFT's

tic;vector_fft=fft(timeseries);tim_est=toc;

vector_fft=[vector_fft,vector_fft];

tim_est = tim_est*ceil((maxfreq - minfreq+1)/freqsamplingrate)   ;

if verbose disp(sprintf('Estimated time is %f',tim_est)),end

% Preallocate the STOutput matrix

st=zeros(ceil((maxfreq - minfreq+1)/freqsamplingrate),n);

% Compute the mean

% Compute S-transform value for 1 ... ceil(n/2+1)-1 frequency points

if verbose disp('Calculating S transform...'),end

if minfreq == 0

   st(1,:) = mean(timeseries)*(1&[1:1:n]);

else

   st(1,:)=ifft(vector_fft(minfreq+1:minfreq+n).*g_window(n,minfreq,factor));

end

%the actual calculation of the ST

% Start loop to increment the frequency point

for banana=freqsamplingrate:freqsamplingrate:(maxfreq-minfreq)

   st(banana/freqsamplingrate+1,:)=ifft(vector_fft(minfreq+banana+1:minfreq+banana+n).*g_window(n,minfreq+banana,factor));

end   % a fruit loop!   aaaaa ha ha ha ha ha ha ha ha ha ha

% End loop to increment the frequency point

if verbose disp('Finished Calculation'),end

%%% end strans function

%------------------------------------------------------------------------

function gauss=g_window(length,freq,factor)

% Function to compute the Gaussion window for 

% function Stransform. g_window is used by function

% Stransform. Programmed by Eric Tittley

%

%-----Inputs Needed--------------------------

%

% length-the length of the Gaussian window

%

% freq-the frequency at which to evaluate

%   the window.

% factor- the window-width factor

%

%-----Outputs Returned--------------------------

%

% gauss-The Gaussian window

%

lambda=0.45;

beta=0;

vector(1,:)=[0:length-1];

vector(2,:)=[-length:-1];

vector=vector.^2;    

vector=vector.*(-factor.*2.*pi^2./freq.^2/lambda.^2-i*2*pi.*vector.*beta./freq);

% Compute the Gaussion window

gauss=sum(exp(vector));

%-----------------------------------------------------------------------

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%

function [ minfreq,maxfreq,samplingrate,freqsamplingrate] =  check_input(minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,timeseries)

% this checks numbers, and replaces them with defaults if invalid

% if the parameters are passed as an array, put them into the appropriate variables

s = size(minfreq);

l = max(s);

if l > 1  

   if verbose disp('Array of inputs accepted.'),end

   temp=minfreq;

   minfreq = temp(1);;

   if l > 1  maxfreq = temp(2);,end;

   if l > 2  samplingrate = temp(3);,end;

   if l > 3  freqsamplingrate = temp(4);,end;

   if l > 4  

      if verbose disp('Ignoring extra input parameters.'),end

   end;

end      

     

   if minfreq < 0 | minfreq > fix(length(timeseries)/2);

      minfreq = 0;

      if verbose disp('Minfreq < 0 or > Nyquist. Setting minfreq = 0.'),end

   end

   if maxfreq > length(timeseries)/2  | maxfreq < 0 

      maxfreq = fix(length(timeseries)/2);

      if verbose disp(sprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d',maxfreq)),end

   end

      if minfreq > maxfreq 

      temporary = minfreq;

      minfreq = maxfreq;

      maxfreq = temporary;

      clear temporary;

      if verbose disp('Swapping maxfreq <=> minfreq.'),end

   end

   if samplingrate <0

      samplingrate = abs(samplingrate);

      if verbose disp('Samplingrate <0. Setting samplingrate to its absolute value.'),end

   end

   if freqsamplingrate < 0   % check 'what if freqsamplingrate > maxfreq - minfreq' case

      freqsamplingrate = abs(freqsamplingrate);

      if verbose disp('Frequency Samplingrate negative, taking absolute value'),end

   end

% bloody odd how you don't end a function

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%

function st_help

   disp(' ')

disp('st()  HELP COMMAND')

disp('st() returns  - 1 or an error message if it fails')

disp('USAGE::    [localspectra,timevector,freqvector] = st(timeseries)')

   disp('NOTE::   The function st() sets default parameters then calls the function strans()')

   disp(' ')  

   disp('You can call strans() directly and pass the following parameters')

   disp(' **** Warning!  These inputs are not checked if strans() is called directly!! ****')

   disp('USAGE::  localspectra = strans(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate,verbose,removeedge,analytic_signal,factor) ')

     

   disp(' ')

   disp('Default parameters (available in st.m)')

disp('VERBOSE          - prints out informational messages throughout the function.')

disp('REMOVEEDGE       - removes the edge with a 5% taper, and takes')

   disp('FACTOR           -  the width factor of the localizing gaussian')

   disp('                    ie, a sinusoid of period 10 seconds has a ')

   disp('                    gaussian window of width factor*10 seconds.')

   disp('                    I usually use factor=1, but sometimes factor = 3')

   disp('                    to get better frequency resolution.')

   disp(' ')

   disp('Default input variables')

   disp('MINFREQ           - the lowest frequency in the ST result(Default=0)')

   disp('MAXFREQ           - the highest frequency in the ST result (Default=nyquist')

   disp('SAMPLINGRATE      - the time interval between successive data points (Default = 1)')

   disp('FREQSAMPLINGRATE  - the number of frequencies between samples in the ST results')

% end of st_help procedure