分数S(Stockwell)变换Matlab代码

离散的FrST

Matlab代码

function FRST= FrST(t,sig,freqlow,freqhigh,alpha,seta)
% Jinshun Shen,Xidian University, March 2022
%jsshen@stu.xidian.edu.cn
% This program can not be used for commercialization without the authorization of its author. 
%FrST 分数S(Stockwell)变换
if size(sig,1)~=1
    sig=sig';
end
if size(t,1)~=1
    t=t';
end

s=sig.*exp(0.5*1i*2*pi*(t.*t)*cot(seta));
st=myst(t,s,freqlow,freqhigh,alpha);
FRST=st.*exp(-0.5*1i*2*pi*(t.*t)*cot(seta));
end

function wcoefs = myst(t,Sig,freqlow,freqhigh,alpha)
% Jinshun Shen,Xidian University, March 2022
%jsshen@stu.xidian.edu.cn
% 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 (nargin < 5)
    alpha=0.05;
elseif (nargin < 4)
    freqhigh=log(length(t));
elseif (nargin < 3)
    freqlow=1;
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=15./fre;
for m = 1:nLevel
% 计算各频率上的ST系数
% 提取频率参数
f= fre(m);
% sigma_f=f*2*pi;
    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] Srivastava H M , Shah F A , Tantary A Y . A family of convolution-based generalized Stockwell transforms[J]. Journal of Pseudo-Differential Operators and Applications, 2020, 11(4).