GST and example.rar_GST 广义S变换_generalized S_s变换广义_广义S变换_时频分析技术；

% 在S变换的基础上，对 gauss=g_window(length,freq,factor) 进行了修改； % 读者可根据需要对其进行其他方式的修改 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, 删除最小二乘拟合抛物线，并在时间序列的边缘放置5％的汉宁锥度。 % This is usually a good idea. %[analytic_signal] if the timeseries is real-valued % 则它将获取分析信号并对其进行ST处理。 %[factor] 局部化高斯的宽度因子，即周期为10秒的正弦波，有一个宽度因子*10秒的高斯窗口。 % 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('没有参数输入'),end st_help t=0;st=-1;f=0; return end % 更改为列向量 (将信号转换为 n*2 的矩阵) if size(timeseries,2) > size(timeseries,1) %矩阵的列数大于行数 timeseries=timeseries'; end % Make sure it is a 1-dimensional array if size(timeseries,2) > 1 error('请输入数据的“向量”，而不是矩阵') return elseif (size(timeseries)==eye(1))==1 error('请输入数据的“向量”，而不是标量') return end % 对输入变量使用默认值 if nargin == 1 minfreq = 0; maxfreq = fix(length(timeseries)/2); %fix:向0靠拢取整 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('输入参数错误：使用默认值'),end minfreq = 0; maxfreq = fix(length(timeseries)/2); samplingrate=1; freqsamplingrate=1; end if verbose disp(fprintf('Minfreq = %d\n',minfreq)) disp(fprintf('Maxfreq = %d\n',maxfreq)) disp(fprintf('采样率（时 域）= %d\n',samplingrate)) disp(fprintf('采样率（频率域）= %d\n',freqsamplingrate)) disp(fprintf('时间序列的长度是： %d points\n',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); % ceil:朝着正无穷方向取整 f = (minfreq + (0:spe_nelements-1)*freqsamplingrate)/(samplingrate*length(timeseries)); if verbose fprintf('The number of frequency voices is %d\n',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 % 如果为0，则绘制幅度谱 if nargout==0 if verbose disp('绘制伪彩色图像'),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, ST Output, 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. % % ST输出的行是频率，列是时间值 % % %----------------------------------------------------------------------- % Compute the length of the data. n=length(timeseries); original = timeseries; if removeedge if verbose disp('用多项式拟合消除趋势'),end ind = [0:n-1]'; r = polyfit(ind,timeseries,2); %%%^^^^^^^^polyfit：线性拟合函数 fit = polyval(r,ind) ; %%%^^^^^^^^polyval：计算多项式 timeseries = timeseries - fit; if verbose disp('用5％汉宁锥度去除边缘'),end sh_len = floor(length(timeseries)/10); %%%^^^^^^^^floor：向负无穷大方向取整 wn = hanning(sh_len); if(sh_len==0) sh_len=length(timeseries); wn = 1&[1:sh_len]; end % 确保wn是列向量，因为时间序列也是列 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 si