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 = FALSE;
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 analy
