广义的S变换
S变换中窗函数是高斯函数 1 2 π σ e − 1 2 σ t 2 \frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{1}{{2\sigma }}{t^2}}} 2πσ1e−2σ1t2,它的形状由方差 σ = 1 f \sigma=\frac{1}{f} σ=f1控制。许多研究表明,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%信号长度大于时间长度,时间补0
t=[t zeros(1,SigLen-TimeLen)];
end
if SigLen < TimeLen%信号长度小于于时间长度,信号补0
Sig=[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.