时频重排技术(Reassignment)原理和Matlab代码

时频重排技术原理

为了提高时频分析工具的分辨率，提出了时频重排技术。它将时频系数重新分配时频重心的位置，以提高时频表示的可读性。重排可以应用于线性的时频分析工具，例如短时傅里叶变换，小波变换等，也可以应用于二次方法，例如魏格纳分布等。

时频表示的能量不够聚集的原因，主要是因为时频系数分布在了信号的瞬时频率的邻域内。那么，只需要将时频系数在时频面上重新分配到估计的瞬时频率处，就可以提高能量聚集性。具体的原理如下图所示，

不同于同步压缩变换，重排是在二维的时频面上重排系数。因此，重排后的时频表示质量，取决于估计的群延迟和瞬时频率。以魏格纳分布为例，它们的估计公式如下：

那么，基于魏格纳分布的重排被定义为

由于重排在二维平面上重分配系数，因此丧失了重构信号的能力。但是，通常重排的时频分辨率更高。

Matlab代码

function [tfr] = RS(x,hlength);
%   Reassignment transform
%	x       : Signal.
%	hlength : Window length.

%	tfr   : Time-Frequency Representation.

%  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.05.13, for my girlfriend.

[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'
%
th=h.*ht;

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

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

va=N/hlength;
    
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));
tfr3(indices,icol)=rSig.*conj(th(Lh+1+tau));
end;

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

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

omega = zeros(round(N/2),tcol);
omega2= zeros(round(N/2),tcol);

for b=1:N
omega1(:,b) = (ft-1)'+real(va*1i*tfr2(ft,b)/2/pi./tfr1(ft,b));
end

for a=1:round(N/2)
omega2(a,:) = t+(hlength-1)*real(tfr3(a,t)./tfr1(a,t));
end

omega1=round(omega1);
omega2=round(omega2);

Ts = zeros(round(N/2),tcol);

% Reassignment step
for b=1:N%time
    for eta=1:round(N/2)%frequency
        if abs(tfr1(eta,b))>0.0001
            k1 = omega1(eta,b);
            k2 = omega2(eta,b);
            if k1>=1 && k1<=round(N/2) && k2>=1 && k2<=N
                Ts(k1,k2) = Ts(k1,k2) + abs(tfr1(eta,b));
            end
        end
    end
end

tfr=Ts/(xrow/2);
end

实验结果

第一幅图是短时傅里叶变换的结果，第二幅图是同步压缩变换的结果，第三幅图是重排的结果。

[1] Auger F , Flandrin P . Improving the readability of time-frequency and time-scale representation by the reassignment method.