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 messageif verbose disp(' '),end% i like a line left blankif nargin == 0if verbose disp('No parameters inputted.'),endst_helpt=0;,st=-1;,f=0;returnend% Change to column vectorifsize(timeseries,2) > size(timeseries,1)timeseries=timeseries';
end% Make sure it is a 1-dimensional arrayifsize(timeseries,2) > 1error('Please enter a *vector* of data, not matrix')returnelseif (size(timeseries)==[11]) == 1error('Please enter a *vector* of data, not a scalar')returnend% use defaults for input variablesif nargin == 1minfreq = 0;maxfreq = fix(length(timeseries)/2);samplingrate=1;freqsamplingrate=1;
elseif nargin==2maxfreq = 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);
elseif verbose disp('Error in input arguments: using defaults'),endminfreq = 0;maxfreq = fix(length(timeseries)/2);samplingrate=1;freqsamplingrate=1;
endif 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 spectrumif nargout==0if verbose disp('Plotting pseudocolor image'),endpcolor(t,f,abs(st))
endreturnfunctionst = 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 removeedgeif verbose disp('Removing trend with polynomial fit'),endind = [0:n-1]';r = polyfit(ind,timeseries,2);fit = polyval(r,ind) ;timeseries = timeseries - fit;if verbose disp('Removing edges with 5% hanning taper'),endsh_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 isifsize(wn,2) > size(wn,1)wn=wn'; endtimeseries(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_signalif 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 pointsif verbose disp('Calculating S transform...'),endif minfreq == 0st(1,:) = mean(timeseries)*(1&[1:1:n]);
elsest(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 pointfor 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 pointif verbose disp('Finished Calculation'),end%%% end strans function%------------------------------------------------------------------------functiongauss=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%vector(1,:)=[0:length-1];
vector(2,:)=[-length:-1];
vector=vector.^2;
vector=vector*(-factor*2*pi^2/freq^2);
% 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 > 1if verbose disp('Array of inputs accepted.'),endtemp=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 > 4if verbose disp('Ignoring extra input parameters.'),endend;endif minfreq < 0 | minfreq > fix(length(timeseries)/2);minfreq = 0;if verbose disp('Minfreq < 0 or > Nyquist. Setting minfreq = 0.'),endendif maxfreq > length(timeseries)/2 | maxfreq < 0 maxfreq = fix(length(timeseries)/2);if verbose disp(sprintf('Maxfreq < 0 or > Nyquist. Setting maxfreq = %d',maxfreq)),endendif minfreq > maxfreq temporary = minfreq;minfreq = maxfreq;maxfreq = temporary;clear temporary;if verbose disp('Swapping maxfreq <=> minfreq.'),endendif samplingrate <0samplingrate = abs(samplingrate);if verbose disp('Samplingrate <0. Setting samplingrate to its absolute value.'),endendif freqsamplingrate < 0% check 'what if freqsamplingrate > maxfreq - minfreq' casefreqsamplingrate = abs(freqsamplingrate);if verbose disp('Frequency Samplingrate negative, taking absolute value'),endend% bloody odd how you don't end a function%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^%functionst_helpdisp(' ')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
