%Author: Yingpin Chen of UESTC and MNNU
%Please cite the folllowing articles if you use the script
%Chen, Y.P., Peng, Z.M.,2012. A novel optimal STFrFT and its application in seismic signal processing.Proceedings of the IEEE International %Conference on ComputationalProblem-Solving, pp. 328-331. %Chen, Y.P., Peng, Z.M., He,Z.H., et al., 2013. The optimal fractional Gabortransform based on the adaptive window function and its %application .Appl. Geophys. 10 (3), 305-313.
function [tfr] = OGTBP(x)
x = x(:);
%============================求最优旋转因子=============================
nfft = length(x);
nOrder=200;
fr = complex(zeros(nfft, nOrder));
for a = 1:nOrder
fr(:,a)=frft(x,2*(a-1)/(nOrder));
end
maxMag = max(abs(fr));
order=0:2/100:2-2/100;
maxMag=maxMag/max(maxMag);
[maxNum, a0] = max(maxMag);
OptOrder=(2*(a0-1)/(nOrder)); %最优阶是OptOrder
sampleInterval=1/sqrt(nfft);
tStart=(-nfft/2)*sampleInterval;
LEN=nfft;
ofr = fr(:, a0);
fofr = (frft(ofr,1)); %对最优阶下的信号做傅里叶变换
normX2 = sum(abs(ofr).^2); %求||x||^2
sampleRate = 1/sampleInterval; %采样率
fRes = sampleRate/LEN; %一格频率代表多少赫兹用fRes表示
tSeries = [tStart:sampleInterval:(tStart+(LEN-1)*sampleInterval)]';
fSeries = [-sampleRate/2: fRes: (-sampleRate/2 + (LEN-1)*fRes)]'; %分数频序号点
itaT = sum((tSeries .* (abs(ofr).^2))/normX2); %中心时间
fLen = length(fSeries);
itaF = sum((fSeries .* (abs(fofr(1:fLen)).^2))/normX2);%中心频率
Tx = (sum(((tSeries - itaT).^2) .* (abs(ofr).^2)))^0.5;%时间标准差
Bx = (sum(((fSeries - itaF).^2) .* (abs(fofr(1:fLen)).^2)))^0.5;%频率标准差
%=======================求TBP最优窗和GTBP最优窗==========================
gama=Bx/Tx;
tStart=(-nfft/2)*sampleInterval;
t = tStart:sampleInterval:tStart + (nfft-1)*sampleInterval;
%求取最优窗
WinTBP__Xp = exp(-pi*gama*(t.^2) );
WinGTBP__X=frft(WinTBP__Xp.',OptOrder).';
tfr = abs(stfrt(1,x, WinGTBP__X,nfft)).^2;
转载本文请联系原作者获取授权,同时请注明本文来自陈颖频科学网博客。 链接地址: https://blog.sciencenet.cn/blog-684084-1058807.html
上一篇:
高级信号处理必看书目 下一篇:
基金申请的几点总结