function [st,t,f] = st(timeseries,minfreq,maxfreq,samplingrate,freqsamplingrate)
% S变换
% Code by huasir @Beijing 2025.1.10
% 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.
%% 函数的输出和输出
% Input
% * timeseries： 待进行S变换的信号向量
% * minfreq： 时频分布结果中的最小频率，对应频率轴的最小值（默认值为0）
% * maxfreq： 时频分布结果中的最大频率，对应频率轴的最大值（默认值为奈奎斯特频率）
% * samplingrate： 两个采样点之间的采样时间间隔（默认为1）
% * freqsamplingrate：频率分辨率（默认为1）
% Output
% * st： Stockwell变换的结果，行对应频率，列对应时刻，
% * t： 包含采样时刻的时间向量
% * f： 频率向量
%% 以下参数可按需调整
% * [verbose]： 如果为真，则打印函数运行中所有的提示信息
% * [removeedge]：如果为真，则删除最小二乘拟合的抛物线，并在时间序列的边缘放置一个5%的hanning锥
% 通常情况下，这是个不错的选择。
% * [analytic_signal]： 如果时间序列是实数值，则取它的解析信号并进行S变换
% * [factor]： 局部化高斯的宽度因子，例如，一个周期为10s的正弦信号具有宽度因子*10s的高斯窗口。
% 通常使用的因子为1，有些情况下为了得到更好的频率分辨率，可以采用3。
% *****All frequencies in (cycles/(time unit))!******
% Copyright (c) by huasir
% $Revision: 1.0 $ $Date: 2025/01/10 $
% 这是保存函数默认值的S变换封装器
TRUE = 1;
FALSE = 0;
%%% 默认参数 [有特殊需求的情况下进行更改]
verbose = TRUE;
removeedge= FALSE;
analytic_signal = FALSE;
factor = 1;
%%% 默认参数设置结果
%% 开始进行输入变量检查
% 首先确保输入的时间序列是有效值，否则，返回帮助信息
if verbose disp(' '),end % i like a line left blank if nargin == 0 % nargin为输入参数的个数，nargin=0，表示无输入if verbose disp('No parameters inputted.'),end st_help t=0;,st=-1;,f=0; return
end % 如果输入数据为行向量的话，将它调整为列向量
if size(timeseries,2) > size(timeseries,1) timeseries=timeseries';
end % 确保输入数据为1维向量，而不是矩阵
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 % 输入变量不全的情况下采用默认值if nargin == 1 %只有1个输入变量minfreq = 0; maxfreq = fix(length(timeseries)/2); samplingrate=1; freqsamplingrate=1;
elseif nargin==2 %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
% 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 > 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
S逆变换函数
function [ts] = inverse_st(st)
% Returns the inverse of the Stockwell Transform of the REAL_VALUED timeseries.
% 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------------------------------------------------
%
% S-Transform matrix created by the st.m function
%-------Optional Inputs ------------------------------------------------
% none
%-------Outputs Returned------------------------------------------------
%
% ts -a REAL-VALUED time series
%--------Additional details----------------------- % sum over time to create the FFT spectrum for the positive frequencies
stspe = sum(st,2); % get st matrix dimensions
[nfreq,ntimes] = size(st);
if rem(ntimes ,2) ~= 0 % odd number of points, so nyquist point is not aliased, so concatenate % the reversed spectrum to create the negative frequencies % drop the DC value negspe = fliplr(stspe(2:nfreq)');
else % even number of points % therefore drop the first point (DC) and the last point (aliased nyqusit freq) negspe = fliplr(stspe(2:nfreq-1)');
end % using symmetry of FFT spectrum of a real signal, recreate the negative frequencies from the positie frequencies fullstspe = [conj(stspe') negspe]; % the time series is the inverse fft of this
ts = ifft(fullstspe);
% and take the real part, the imaginary part will be zero.
ts = real(ts);
