广义的S变换

广义的S变换

S变换中窗函数是高斯函数$\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{1}{2\sigma }{t}^2}$，它的形状由方差$\sigma =\frac{1}{f}$控制。许多研究表明，S变换中窗函数的宽度并不是最优的，导致S变换对非平稳信号刻画地不够细致。因此，通过引入额外的参数来控制高斯窗的形状，提出了广义的S变换（自适应的S变换），提升S变换的时频分辨率。本质上，广义的S变换是通过改变高斯函数的方差来实现的。对比广义的S变换和S变换代码，很容易可以看出这一结论。

广义的S变换代码

function wcoefs = G_myst(t,Sig,freqlow,freqhigh,alpha,lamda,p,q)
%==================%
% Stockwell transform ver 1.0 (2021-09-24)
% Copyright (c) by Jinshun Shen
%This program can not be used for commercialization without the authorization of its author
%======输入======%
%t:时间序列
% Sig: 输入信号
%freqlow,freqhigh:频率范围
%alpha:频率分辨率
%======输出======%
% wcoefs: ST变换计算结果if (nargin <= 2)
error('At least 2 parameters required');
end
if size(t,1)>1%t是列向量,则转置t=t';
end
if size(Sig,1)>1%Sig是列向量,则转置Sig=Sig';
end
% 信号的长度
SigLen = length(Sig);
% 时间的长度
TimeLen = length(t);
dt=t(2)-t(1);
fre=freqlow:alpha:freqhigh;%产生频率范围
% time=1:TimeLen;
if SigLen > TimeLen%信号长度大于时间长度，时间补0t=[t  zeros(1,SigLen-TimeLen)];
end
if SigLen < TimeLen%信号长度小于于时间长度，信号补0Sig=[Sig  zeros(1,TimeLen-SigLen)];
end
% 总频率数量
nLevel=length(fre);
% 分配计算结果的存储单元
wcoefs = zeros(nLevel,TimeLen); 
temp=zeros(1,TimeLen);
% sigma_f=lamda*((2*pi*fre).^p)+q;
sigma_f=lamda*((fre).^p)+q;
for m = 1:nLevel
% 计算各频率上的ST系数
% 提取频率参数
f= fre(m);for n=1:TimeLen%计算高斯窗函数Gauss_st=(1/(sqrt(2*pi)*sigma_f(m)))*exp(-0.5*((n*dt-t).^2/(sigma_f(m))^2)).*exp(-1i*2*pi*f*t);temp(n)=trapz(t,Sig.*Gauss_st);%效率低，精度稍高end
wcoefs(m,:)=temp;
end
end

参考文献

[1] R. G. Stockwell, L. Mansinha, and R. Lowe, “Localization of the
complex spectrum: the S transform,” IEEE transactions on signal
processing, vol. 44, no. 4, pp. 998–1001, 1996.

[2] Wang Q, Gao J, Liu N and Jiang X. High-resolution seismic time-frequency analysis using the synchrosqueezing generalized S-transform [J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15 (3): 374-378 (SCI: 000426272400012).

[3] Wang L , Meng X . An adaptive Generalized S-transform for instantaneous frequency estimation[J]. Signal Processing, 2011, 91(8):1876-1886.