S变换matlab实现

S变换函数

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);