【时频分析】广义线性chirplet变换【附MATLAB代码】

文章来源：微信公众号：EW Frontier

1.摘要

时频分析方法是一种刻画信号时变特征的有效工具，在相当长的一段时间内受到了广泛的关注。随着TF算法的发展，许多先进的方法被提出，可以提供更精确的TF结果。但是，不可避免地会引入一些限制。本文介绍了一种新的TFA方法，称为广义线性线调频小波变换（GLCT），它可以克服现有TFA方法存在的一些局限性。在数值和实验验证中，通过与现有的TFA方法的比较，GLCT的一些优点被证明，包括良好的表征具有明显的非线性特征的多组分的信号，独立于数学模型和初始TFA方法，允许感兴趣的组分的重建，和对噪声不敏感。

2.广义线性chirplet变换（GLCT）

2.1.GLCT的动机

首先，我们用一个数值例子来说明我们所提出的方法的动机。该模拟信号借用自参考文献[22]，如下所示：

其以100 Hz的采样频率采样，并且其IF是时变的。图1（a）和（B）中所示的TF表示分别由STFT和LCT（啁啾率c =4π）生成。从本质上讲，STFT可以看作是线性调频为零的LCT。从图1（a）和（B）中的TF结果可以看出，当LCT的啁啾率接近信号的真实IF时，TF表示将具有高能量集中。相反，如果LCT的啁啾率远离信号的真实IF，则TF表示将严重拖尾。如上所述，当所考虑的信号包含时变频率分量时，任何单独的LCT都不可能生成具有令人满意的能量集中的TF表示。因此，有必要开发替代方法来解决这一问题。一个直观的解决方案的动机是通过结合的一部分，高能量集中的结果中产生的几个LCT与不同的啁啾率。并且新的TF表示将具有比任何单独的LCT更好的性能，如图1（c）所示。因此，出现了两个需要解决的问题。第一个问题是如何确定LCT的完整啁啾率，以确保它们能够尽可能精确地描述信号的时变频率。第二个问题是如何确定标准，以选择合适的TF表示由不同的LCT产生。在接下来的章节中，我们将给予解决这两个问题的详细理论分析。

图1.通过（a）STFT、（b）LCT和（c）GLCT表示TF。

2.2.GLCT理论

根据Ville的理论，考虑时变IF的解析信号可以写为

其中A（t）是幅度，其IF可以写为

在短时间μ内，信号s（t）的中频可近似视为一个线性方程，可采用泰勒展开式进行展开

其中忽略提醒项，φ （t′）是时间t′处的IF，φ ′（t ′）是时间t′处的IF的一阶导数。

给定一个窗口w （u−t′）来截断信号，则IF φ （t′）中w （u−t′）*s( u)的FT可以计算为：

当量(5)表示由于调制元件e的存在，w （u−t′）*s( u)在其IF中的FT幅度低于谐波分量。实际上，Eq.(5)是分析信号在其IF中的STFT幅度。这就是为什么调制部分的TF表示能量在中频周围大面积扩散，而谐波部分的TF表示能量可以集中在幅度较高的中频。

给定标准STFT公式，

为了消除调制分量的影响，需要引入解调算子。考虑到线性调制比φ  ′（t′）是时间的函数，即时变的，解调算子也应该是时变的。然后，考虑时变解调算子的信号s （t）的STFT可以写为：

对于Eq.(7)，IF φ（t′）中s （t）的STFT幅度可通过下式计算：

从等式(8)可以看出，如果解调算子与分析信号的调制成分一致，则中频φ t（′）中的STFT幅度将达到最大值，此时TF表示将集中在能量集中度高的中频。然而，在大多数实际情况下，信号的中频特性并不能预先知道，这就导致了解调后的算子e很难准确确定，尤其是当信号中存在多分量时。在本文中，我们不打算识别信号的中频特征，并提出了一种替代和易于编程的方法。

一个可行的解决方案是使用一系列离散解调算子来近似最佳解调算子e。给出考虑离散解调算子的STFT公式，如下

对于每个TF点（t′，ω），如果离散解调算子e接近信号的调制元素，则其IF周围的TF表示将具有高能量集中，并且幅度丨S（t′，ω，c）丨将达到所有值中的最大值。然后，对于每个TF点，根据丨S（t′，ω，c）丨的幅度，我们可以获得最佳幅角c，

并且我们提出的TFA方法的TF表示可以获得为

频谱可以定义为：

下一个需要解决的问题是如何确定离散解调算子e。事实上，Eq。(9)与线性线调频小波变换等价，这就是我们提出的方法被称为广义线性线调频小波变换（GLCT）的原因。根据参考文献的研究结果。[22，25]，通过引入e ，会对TF结果产生旋转效应，旋转度为arctan（−c ）。对于信号，如果采样时间（Ts）和采样频率（Fs）被确定，则它们将确定TF平面t∈（0，Ts）, f∈（0，Fs/2）。因此，我们可以引入一个参数α，它可以在这个TF平面中旋转，

这使得解调算子e能够描述信号中所有可行的调制元素，则等式(9)可以重写为

假设参数α有N个值，可以将TF平面平均分成N+1个截面，

从上面的等式可以看出，我们提出的方法只比STFT多引入一个参数，即N。如果N =1，GLCT将退化为STFT。

我们用一个图表进一步说明了GLCT（N = 3）的原理。在图2中，实线表示真实IF，虚线框表示分析窗口，灰色区域表示TF表示的能量分布。在图2（a）和（B）中，由于使用了单独的啁啾窗口，如果解调比接近信号的IF，则TF表示的能量将集中在具有高TF分辨率的IF上。然而，如果解调比与信号的IF显著不同，则TF表示的能量将在大的区域中扩散。在图2（c）中，通过在几个LCT产生的结果中选择具有高能量集中的部分，GLCT可以在IF周围实现比任何单个LCT更高的能量集中，因此它可以为我们提供更精确的TF结果。

图2.（a）STFT、（b）LCT和（c）GLCT的TF能量分布。

3.数值验证

在这一节中，我们使用了几个数值例子，其中包括一个单分量信号和两个多分量信号来说明GLCT的性能与其他TFA方法相比。

3.1.单分量信号

第一个信号被认为是

通过不同TFA方法获得的结果如图3所示。如图3（a）所示，由WT生成的TF表示的能量在IF周围扩散，这导致差的TF分辨率。在图3（B）中，由于意外交叉项存在，由WVD产生的TF表示不能给予关于信号的有用信息。对于图3（c）中所示的STFT的结果，TF表示的能量在谐波部分中比在调制部分中更高。可以看出，由于Heisenberg不确定性的限制，WT和STFT都不能同时在时域和频域上实现高TF分辨率，特别是在分析调频信号时。

图三.通过（a）WT、（b）WVD、（c）STFT、（d）AD、（e）SCT、（f）GLCT（N1/4 5）、（g）GLCT（N1/4 7）和（h）GLCT（N1/4 15）的TF表示。

原子分解（AD）方法可以实现稀疏TF表示[30]，但高度依赖于所选原子字典[31]。对于数值信号，选择“Symmlet”作为原子函数，TF结果如图3（d）所示。可以看出，对于非平稳信号，所选原子不是一个好的选择。样条线调频小波变换（SCT）是GPTFA [23]的一种特殊情况，用于分析该信号。首先，需要在原始短时傅里叶变换的基础上识别中频特征，然后在TF平面上构造数学模型来表征信号，如图3（e）所示。时变TF特性基本得到表征，但需要额外的过程来估计IF。

对于图3（f）、（g）和（h），它们是通过GLCT以不同参数生成的结果，即分别为N= 5、7和15，并且窗口长度（WL）被设置为100个点。结果表明，随着参数N的增加，GLCT产生的TF表示具有更高的能量集中度，更接近理想的TF表示。

实际上，当N= 7和15时，GLCT产生的TF表示对于表征信号的真实IF特征是令人满意的。

3.2.中频特征估计

信号的中频是一个重要的特征，它可以为我们提供关于分析对象的重要信息。在这一小节中，我们测试了GLCT检测信号中频特征的能力，以验证GLCT的抗噪声性能。并选择短时傅里叶变换作为比较基准，在方程（16）的信号中加入不同信噪比的噪声。而中频可以通过TF表示的峰值数据来检测。检测到的误差可以通过下式计算：

其中，IF  e（t）表示估计的IF，IF （t）是真实的IF。图4示出了通过两种TFA方法估计IF的计算误差。可以看出，通过GLCT（N = 7和WL = 100）估计的IF比通过STFT估计的IF对噪声更鲁棒。为了更清楚地说明检测过程，我们列出了SNR =10 dB和SNR =0 dB下的TF表示和估计IF。

图4.STFT和GLCT检测中频的误差。

如图5（a）和图5（B）所示，当加入SNR= 10 dB的噪声时，STFT检测到的IF在真实IF附近，并且估计误差很低。然而，当加入信噪比为0 dB的噪声时，STFT检测到的中频与真实中频发生了较大的偏差，特别是在调制部分，如图6（a）和（b）所示。这是由于信号中调制部分的幅值较低，比谐波部分更容易受到噪声的影响。GLCT的检测结果示于图1A和1B中。如图5（c）和（d）以及图6（c）和（d）所示，当添加SNR的=10 dB和SNR=0 dB的噪声时，估计的IF都显示出与真实IF的令人满意的精度。

3.3.信号重构

TFA的主要目标可以表述为信号中感兴趣分量的识别和重建。因此，重建的量化是评价TFA方法的另一个重要指标。对于传统的方法，如STFT和WT，有两种可行的重建类型，直接重建和脊重建。直接法的重构是通过选择TF表示中的感兴趣区域，利用TF表示的逆公式提取指定分量，是一种在大多数情况下应用广泛的方法。然而，直接重建方法不适合GLCT。脊重建已经在许多文献中进行了研究，这表明可以使用脊点处的TF表示数据来恢复信号[26，27]。

图8.用STFT、WT和GLCT重建5- 6.5s的信号。(For在这幅图中对颜色的引用的解释，读者可以参考本文的网络版本。

对于GLCT的TF表示，当感兴趣分量的脊被估计时，信号可以通过重建

其中R（.）表示取真实的部分，G S（t ，ω）是GLCT的TF表示，并且r（t）是检测到的脊。

在这一小节中，我们测试的能力GLCT重建的信号方程。(16)在不同的噪音水平下对于GLCT，我们使用从峰值数据检测到的IF作为岭参数r（t）。我们选择短时傅立叶变换和小波变换作为比较基准，利用直接法重构信号，并以重构信号的信噪比作为评价指标。评价结果如图7所示。可以看出，GLCT在低噪声水平下具有最好的重建能力。对于噪声信号（SNR = 0 dB），我们列出了5-6.5 s的重建结果（参见图8，蓝色实线为原始信号，红色虚线表示重建信号），这表明GLCT在所选TFA方法中提供了最佳重建。这是由于附加噪声在TF平面上的传播范围很广，所以TF平面上的重构区域越小，引入到重构信号中的噪声就越少，GLCT的公式只考虑了IF上的重构区域，因此它应该具有最好的重构性能。

3.4.多分量信号

具有多分量的第二数值信号被认为是

其IF绘制在图9（a）中。可以看出，分量s1（t）、s2（t）具有相同的调制分量但具有不同的初始频率，并且s3（t）具有与其它分量不同的调制分量。为了与其他方法进行比较，我们列出了由不同方法生成的TF表示，例如WT，STFT，GPTFA，SST和GLCT（N =7和WL = 100）。

图9.(a)通过（b）WT、（c）STFT、（d）GPTFA、（e）SST和（f）GLCT表示真IF、TF。

通过不同方法生成的TF表示，例如WT、STFT、GPTFA、SST和GLCT（N = 7和WL = 100）。传统的线性时域分析方法，短时傅里叶变换和连续小波变换，都是建立在所考虑的信号是短时分段平稳的假设上的。如图9（B）和（c）所示，当信号的IF恒定时，通过WT和STFT产生的TF表示可以集中在具有高能量集中的真实IF上。然而，对于调制部分，TF表示的能量在IF周围的大区域中扩散，这导致低分辨率。因此，传统的方法只能对时变中频信号提供精确的分析结果。

对于GPTFA，通过识别信号的中频特征并构造非线性解调算子，TF表示可以对中频特征被良好识别的信号实现高TF分辨率。然而，对于在中频特征明显的分量情况下，由于所选参数与分量之间的不匹配，GPTFA甚至可能提供比STFT更差的TF表示。对于图9（d）所示的TF结果，选择s1（t）、s2（t）的数学模型作为GPTFA的参数，因此可以清楚地表征这两个信号。但是，s3（t）的TF结果显示出较差的能量集中。因此，GPTFA仅限于分析具有不同中频特征的多分量信号。因此，GPTFA的主要任务集中在将多分量信号分离为多个单分量信号，并构建精确的数学模型来识别单分量信号的IF特征[28，29]。

SST作为一种后处理TFA方法，必须限制初始TFA方法的缺点。存在两种可行的SST解决方案，基于WT的SST和基于STFT的SST。图9（e）中的结果由后一种方法生成。可以看出，谐波部分在所有TFA方法中具有最高的能量集中，然而，调制部分拖尾严重。因此，为了提高SST的性能，一个可行的解决方案是改进初始TFA方法的结果[14]。

对于图9（f）中的结果，GLCT提供了令人满意的TF表示，其清楚地表征了所有分量。能量分布均匀地集中在真中频上。

3.5.具有交叉中频的多分量信号

在3.4节中，数值多重分量具有良好分离的IF。但是，在实际应用中经常会遇到中频交叉的情况，需要考虑。给出三个数字源，

图10.（a）S12和（B）S13的波形，（c）S12和（d）S13的TF表示。

实际上，分量s 1（t）是等式（1）的信号。分量s 2（t）和s 3（t）是纯谐波信号，其具有相同的频率但不同的相位。两个多分量信号分别被认为是s12=s1+s2和s13=s1+s3。因此，在TF平面中，信号s12和s13应该在（4.3 s和25 Hz）处具有相同的交叉TF点。为了更清楚地展示GLCT的性能，我们列出了s12和s13的波形，如图10（a）和（B）所示。可以看出，在交叉时间4.3s处，s13的幅度低于s12的幅度。这是由于s1和s3的叠加对s13的幅度产生了减小的影响。而s1和s2的叠加具有增加的效果。除了交叉的TF点外，其他点具有相同的幅度。根据GLCT的TF结果（见图10（c）和（d））（N = 7和WL = 100），它们清楚地反映了这一现象。

4.参数选择

对于GLCT，如何选择合适的参数，如N，窗口类型，窗口长度，是产生良好的TF结果的关键。因此，有必要分析不同参数对TF表示的影响。

4.1.参数N的选择

从图3中的数值结果可以看出，在选择较大N的情况下，TF表示在IF中具有较高的集中能量。因此，可以更准确地实现IF的估计。为了清楚说明参数N的影响，我们使用中频估计的误差来评估GLCT的性能。Eq(16)的信号。在不同的噪声水平（10 dB，5 dB和0 dB）的情况下。估计结果见图11。可以看出，随着N的增加，中频估计的误差减小。对于N＞11，误差减少不明显。然而，计算成本将增加。

4.2.窗型选择

选择合适的窗口将有利于TF结果，因此需要确定最佳的最优窗口。常用的窗有汉明窗、汉宁窗、高斯窗和矩形窗等。(16)在不同的噪声水平下，并评估GLCT的性能。选择峰值数据与IF之间的误差作为评价指标，如图12所示。可以看出，高斯窗在高信噪比水平下具有最好的抗噪声鲁棒性，但在低信噪比水平下表现出较差的性能。汉宁窗在低信噪比条件下具有较好的中频估计效果。因此，高斯窗更适合于高信噪比的信号，而汉宁窗更适合于低信噪比的信号。

图12.采用不同窗函数检测中频的误差。

4.3.窗口长度的选择

由于GLCT建立在线性TFA的基础上，因此它受到时间分辨率和频率分辨率之间不确定性关系的限制。为了获得高的频率分辨率，我们应该使用长窗口。根据等式(4)在短时间窗内，采用一阶泰勒展开式对中频进行分析，表明中频在短时间内近似为线性。然而，在这短时间内，如果IF是强非线性的，则泰勒余项将很大，这可能导致这种近似的大误差。因此，对于强非线性中频，我们需要选择一个短的窗口，以减少近似误差。

我们利用一个数值例子来说明窗口长度的TF表示的效果。我们把Eq的信号(16)分析并列出不同窗长（WL = 50、120、250、400点，N = 7）GLCT的TF结果，如图13所示。可以看出，当窗口长度为50点时，TF表示（参见图13（a））具有低频率分辨率。当窗口长度增加时，频率分辨率提高。但是对于400点的窗口长度，由于泰勒展开的近似误差较大，TF结果不能提供令人满意的TF结果（参见图13（d））。因此，建议在频率分辨率和窗口长度之间进行权衡。

图13.TF表示使用不同的窗口长度（a）50、（B）120、（c）250和（d）400。

5.MATLAB代码及仿真结果
 

%A numerical signal to validate the ability of characterizing the signal with time-varying IF and reconstructing it.
%This signal is noise-free, so you can add noise with different level to test its ability.
SampFreq = 100;
hlength=100;
​
t = 0 : 1/SampFreq : 14-1/SampFreq;
Sig = [sin(2*pi*(25*t + 10*sin(t))) ];
​
time=(1:1400)/100;
fre=1/700:50/700:50;
​
%N=7.
N=7;
tfr = GLCT((Sig)',N,SampFreq,hlength);
figure
imagesc(time,fre,abs(tfr).^2);
axis xy
ylabel('Freq / Hz');
xlabel('Time / Sec')
title('GLCT,N=7');
%N=15.
​
N=15;
tfr = GLCT((Sig)',15,100,100);
figure
imagesc(time,fre,abs(tfr).^2);
axis xy
ylabel('Freq / Hz');
xlabel('Time / Sec')
title('GLCT,N=15');
​
%Reconstruction of Signal
[m, n] = size(Sig);
[v, I] = max(abs(tfr),[],1);
for j=1:n
rsig(j)=real(tfr(I(j),j));
end
​
hlength=hlength+1-rem(hlength,2);
h = tftb_window(hlength);
h=h/norm(h);
​
rsig=rsig/(sum(h)/2);
%'rsig' is the reconstructed signal form GLCT, by ridge reconstruction method.
%'cor' is the correlation coffecient between original signal and recovered signal.
cor=corr(Sig',rsig');
%plot the original signal and recovered signal.
figure
plot(Sig);hold on;
plot(rsig,'r');
title('Original signal (blue). Recovered signal (red).');
​
​

function [tfr,tt,ff] = GLCT(x,N,fs,hlength);
%      General Linear Chirplet Transform.
%  x      : Signal.
%  N      : Number of calculated LCT with different chirp rate.
%  fs     : Sample Frequency .
%  hlength: Length of window function.
​
%  tfr    : Time-Frequency Representation.
%  tt     : Time.
%      ff     : Frequency.
%
%  This program is free software; you can redistribute it and/or modify
%  it according to your requirement.
%
%  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 Gang Yu in Shandong University at 2015.4.28.
​
​
​
​
[xrow,xcol] = size(x);
if (nargin < 3),
error('At least 3 parameter is required');
end;
Siglength=xrow;
​
if (nargin < 4),
hlength=floor(Siglength/4);
end;
​
hlength=hlength+1-rem(hlength,2);
h = tftb_window(hlength);
​
%t=1:xrow;
%[trow,tcol] = size(t);
​
[hrow,hcol]=size(h); Lh=(hrow-1)/2; 
​
h=h/norm(h);
​
Frelength=round(Siglength/2);
​
slope=(-pi/2+pi/(N+1)):pi/(N+1):(pi/2-pi/(N+1));
Spec=zeros(N,Frelength,Siglength);%%%
​
t=Siglength/fs;
​
for i=1:N
[Spec(i,:,:)] = LCT(x,tan(slope(i))*(fs/2)/t,fs,h);
end 
​
GLCTspec=zeros(Frelength,Siglength);
​
for i=1:Frelength
for j=1:Siglength
absSpec=abs(Spec(:,i,j));
index=find(max(absSpec)==absSpec);
GLCTspec(i,j)=Spec(index(1),i,j);
end
end
​
tfr=GLCTspec;
​
if (nargout==2),
tt=(1:Siglength)/fs;
elseif (nargout==3),
tt=(1:Siglength)/fs;
ff=0:fs/Siglength:(fs/2)-(fs/2)/Siglength;
end;
​