同步提取变换原理和Matlab代码

同步提取变换(SET)

同步提取变换（synchroextracting transform, SET），它属于短时傅立叶变换（short-time Fourier transform, STFT）的后处理过程。与经典的时频分析(time-frequency analysis, TFA)方法相比，SET可以生成能量更集中的时频表示(time-frequency representation, TFR)，并允许信号重建。SET受到了同步压缩变换(synchrosqueezing transform , SST)和理想的时频表示理论的启发。

SST是将系数压缩到估计的瞬时频率处，而SET只保留与信号的时变特征最相关的STFT系数，并去除最模糊的TF能量，从而可以大大提高新的TF表示的能量浓度。换句话说，SET丢弃了部分STFT系数，以获得高能量聚集性的TFR。它的缺点也是显然的，损失了重构的精度。以下图为例，SET提取${\omega }_0$附近处的系数，它的瞬时频率估计方法与SST类似。

Matlab代码

function [IF Te tfr] = SET(x,hlength);
%   Synchroextracting Transform
%	x       : Signal.
%	hlength : Window length.

%   IF   : Synchroextracting operator representation.
%   Te   : SET result.
%   tfr  : STFT result
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
%
%   Written by YuGang in Shandong University at 2016.5.13.

[xrow,xcol] = size(x);

N=xrow;

if (xcol~=1),
 error('X must be column vector');
end;

if (nargin < 2),
 hlength=round(xrow/8);
end;

t=1:N;
ft = 1:round(N/2);

[trow,tcol] = size(t);

hlength=hlength+1-rem(hlength,2);
ht = linspace(-0.5,0.5,hlength);ht=ht';

% Gaussian window
h = exp(-pi/0.32^2*ht.^2);
% derivative of window
dh = -2*pi/0.32^2*ht .* h; % g'

[hrow,hcol]=size(h); Lh=(hrow-1)/2;

tfr1= zeros (N,tcol);
tfr2= zeros (N,tcol);

for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
rSig = x(ti+tau,1);
tfr1(indices,icol)=rSig.*conj(h(Lh+1+tau));
tfr2(indices,icol)=rSig.*conj(dh(Lh+1+tau));
end;

tfr1=fft(tfr1);
tfr2=fft(tfr2);

tfr1=tfr1(1:round(N/2),:);
tfr2=tfr2(1:round(N/2),:);

va=N/hlength;
IF=zeros(round(N/2),tcol);
tfr=zeros(round(N/2),tcol);
E=mean(abs(x));

for i=1:round(N/2)%frequency
for j=1:N%time
     if abs(tfr1(i,j))>0.8*E%if you are interested in weak signals, you can delete this line.
         %if abs(1-real(va*1i*tfr2(i,j)/2/pi./tfr1(i,j)))<0.5
         if abs(-real(va*1i*tfr2(i,j)/2/pi./tfr1(i,j)))<0.5
         IF(i,j)=1;
         end
     end
end
end
tfr=tfr1/(sum(h)/2);%the amplitude of tfr result has been pre-rectified.
Te=tfr.*IF;
end

实验结果

第一幅图是STFT结果，第二幅图是SET的结果。

[1] Gang Y , Yu M , Xu C . Synchroextracting Transform[J]. IEEE Transactions on Industrial Electronics, 2017, 64(10):8042-8054.