【时频变换】多项式Chirplet变换（PCT）及其在瞬时频率估计中的应用【附MATLAB代码】

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

文章信息：

Polynomial Chirplet Transform With Application to Instantaneous Frequency Estimation

Digital Object Identifier 10.1109/TIM.2011.2124770

1.摘要

在本文中，一个新的时频分析方法称为多项式线性调频波变换（PCT）是通过扩展传统的线性调频波变换（CT）。通过在CT中使用多项式函数代替线性啁啾核，PCT可以产生具有连续瞬时频率（IF）的宽范围信号的具有优异集中度的时频分布。在此基础上，提出了一种有效的中频估计算法，并将其应用于含非线性啁啾、受高斯噪声严重污染的信号和转子试验台振动信号的中频估计，验证了该算法的有效性。

索引项-线性调频小波变换（CT），瞬时频率（IF），多项式线性调频小波变换（PCT），短时傅立叶变换（STFT），时频分析。

2.引言

瞬时频率（IF）的概念是由卡森和弗莱[1]于1937年提出的，用于分析单分量调频信号，其后由货车德尔波尔[2]和伽柏[3]进一步发展。现在，Ville [4]提出了广泛认可的IF定义，统一了卡森和Fry [1]和Gabor [3]的工作。根据Ville的定义，信号的IF可以从其解析信号的相位的导数计算。

IF概念在研究各种信号中发挥了重要作用，从雷达[5]，声纳[6]，生物医学工程[7]，地震调查[8]，语音和音乐[9]到汽车信号[10]，IF通常用于表征信号的重要物理参数。已经开发了各种IF估计方法;Boashash [11]，[12]对这些进行了全面的概述，其中介绍了IF概念本身的解释以及1992年之前发表的关于IF估计方法的文献。一般来说，中频估计方法可以分为五类：1）基于相位差的中频估计器; 2）过零中频估计器; 3）基于线性预测滤波器的自适应中频估计器; 4）基于时频分布矩的中频估计器;和5）基于时频分布峰值的中频估计器。在五类IF估计器中，使用TFD可以产生对噪声更可靠和更鲁棒的结果。因此，基于时域有限差分的方法受到了广泛的关注，并得到了较快的发展。

基于TFD的方法的能力取决于TFD在时频平面中将信号的能量集中在IF处及其周围的特性。三种方法，即短时傅立叶变换（STFT）[13]，连续小波变换（CWT）[14]，[15]和Wigner-Ville分布（WVD）[16]，通常用于产生信号的TFD。短时傅里叶变换和连续小波变换本质上是一种线性变换，其特征在于在时频平面上的静态分辨率，其被细分为恒定面积的基本单元。然而，由于Heisenberg-Gabor不等式的限制，STFT和CWT都不能在时域和频域都获得很好的分辨率，好的时间分辨率必然意味着差的频率分辨率。因此，这两种方法得到的时频表示只有在一定的时频分辨率下才有意义。因此，STFT和CWT永远不能产生信号的真实时间-频率模式，因此只能为信号的IF提供有限测量精度的估计，特别是当IF轨迹是时间的非线性函数时。为了提高估计精度，已经提出了重新分配技术[17]。WVD是一种二次变换，可以实现对无噪声信号的高精度估计。当信号具有线性频率规律和恒定幅度时，其WVD将沿着线性IF轨迹减少到一行δ函数沿着。然而，在噪声信号或非线性IF规律的情况下，WVD基于峰值的估计器将是有偏的，并且已经表明，对于估计的IF结果，偏差与方差的折衷是不可避免的。由IF非线性引起的偏置与滞后窗口长度的幂成比例，而由噪声引起的方差是滞后窗口长度的递减函数。已经做出努力来提高WVD的能力，以实现对非线性IF的更好的估计，主要是通过调整窗口长度。Stankovic和Katkovnik [18]提出了一种方法，WVD具有自适应窗口长度，并且同一作者还提出了一种通过使用具有变化和数据驱动窗口长度的WVD的方法[19]。基于TFD峰值的算法也应用于其他WVD表示，包括多项式WVD [20]，[21]，其被证明对于非线性IF轨迹是无偏的，以及伪和平滑伪WVD [22]，其旨在抑制由时频平面中的噪声引起的交叉项。

除了前面提到的TFD分析方法外，线调频小波变换（CT）[23]是另一种TF方法，专门用于分析具有线性IF规律的线性调频信号。当所考虑的信号的IF轨迹是时间的非线性函数时，CT不再能够保证分析精度的提高，因此Chassande-Mottin和派[24]建议使用线性调频波链来处理本质上是准周期FM信号的引力波信号，Angrisani和D 'Arco [25]提出了一种CT的修改版本，通过将额外的曲率参数引入到传统CT中，这被证明对于具有以较高动态为特征的IF的信号特别有效。本文提出了一种新的CT称为多项式CT（PCT）。基于PCT，然后开发的程序来估计具有高度非线性的IF轨迹的信号的IF。PCT是通过用多项式非线性IF法的新核取代CT中线性IF法的线调频核而开发的。在数学中，维尔斯特拉斯逼近定理[26]保证了有界区间上的任何连续函数都可以用多项式以任何精度一致逼近。因此，PCT能够对各种信号进行高精度分析，其中IF轨迹是时间的任何连续函数。

本文的组织结构如下。在简要描述了传统CT的理论基础之后，第二节给出了关于PCT的细节。基于PCT的IF估计过程在第三节中建立。在第四节中，新提出的中频估计器的有效性进行了验证，通过将其应用到几个信号。结论见第五节。

3.理论分析

A.常规CT

图. 1.传统CT的图示（--目标啁啾信号的IF定律;--频率旋转后;--频移后）。

信号s（t）∈ L2（R）的CT [23]，[27]定义为：

其中z（t）是由希尔伯特变换H生成的s（t）的分析信号，即，z（t）= s（t）+ jH[s（t）]，并且φ（t0，α，σ）（t）是由下式给出的复窗口：

这里，参数t0，α ∈ R分别代表时间和啁啾率; w ∈ L2（R）表示非负，和归一化的真实的窗口，通常取为高斯函数，表示为

根据（1）给出的CT的定义，CT可以被解释为分析信号的STFT乘以复窗口φ（t0，α，σ）（t）。CT的定义也可以写为

显然，ΦR α（t）是一个频率旋转算子，它在时频平面中将分析信号z（t）旋转一个角度θ，其中tg（θ）= −α; ΦM α（t，t0）是一个频移算子，它将ω处的频率分量重新定位到ω + αt0; A（t0）是一个具有模的复数|A（t0）|=1。在时频分析中，它是TFD的模|CT（t0，ω，α; σ）|这通常是有意义的，因此，CT的定义可以简化为

从这个定义可以看出，CT可以分解为一系列算子：1）在时频平面上将所考虑的信号旋转一个角度arctan（−α）; 2）将信号移动一个频率增量αt0; 3）使用窗口w（σ）进行STFT。这个过程可以由图1来说明，其中对象啁啾信号具有ω0 + λ 0 t的IF定律。

图2.CT生成的TFD（窗口大小= 512）。(a)λ0 =0。(b)λ0 =5π。

CT的频率分辨率可以很容易地从图1中确定。如果啁啾率参数α =0，CT等于标准STFT，由具有时间带宽σ的高斯窗选择的信号部分的频率带宽为σλ0，则在时频平面中，该信号部分将被表示为频率带宽σλ0 + 1/σ（1/σ是高斯窗的频率带宽）和时间带宽σ的单元。在α =0的情况下，所选信号部分的频率带宽将减小到σ| λ0 − α| +1/σ，时间带宽保持不变。可以看出，当α = λ0时，频率带宽将具有最小值1/σ;这意味着CT产生的TFD具有最佳浓度。因此，我们认为，|CT（t0，ω，α; σ）|在（ω，α）=（ω0，λ0）处具有全局最大值。

考虑一个例子，其中CT被应用于由两个分量组成的信号，其中一个分量的IF定律为10 + 2.5t（Hz），另一个分量的IF定律为12 + 2.5t（Hz），即，

以200 Hz的采样频率对信号进行采样。使用两种不同的啁啾速率，即，λ0 =0且λ0 =5π。结果分别示于图2（a）和（B）中。可以看出，当λ0 =0时，CT退化为正常的STFT，并且由于大的频率带宽，在图2（a）所示的TFD中不能分离两个分量;相反，如图2（b）所示，这两个分量在这种情况下在时间-频率平面中清楚地显示出来，并且TFD的集中度显著提高，因为CT中使用的啁啾速率参数被适当地选择以匹配所考虑的信号的啁啾速率。

B.PCT

如第II-A节中的例子和其他研究人员提出的一些应用所示，当适当选择啁啾率时，常规CT将使中频轨迹是时间的线性函数的信号的TFD具有良好的集中性。然而，当信号的IF轨迹不完全是时间的线性函数时，则常规CT不可能跟踪该演变，与图2（B）中一样接近的信号IF的时间。这是特别真实的，在存在的IF特征在于高度非线性动力学的线性近似是不足以局部表示的IF轨迹。为了提高传统CT在分析具有非线性中频轨迹的信号时的效率，本文提出了一种称为PCT的改进型CT。PCT的定义如下：

图3.图为PCT。

因为维尔斯特拉斯逼近定理[25]保证了在封闭和有界区间上的任何连续函数都可以在该区间上被多项式一致逼近到任何精度，例如，通过使用最小二乘法，图4中所示的阶梯函数可以很好地被15阶多项式函数逼近，使得

图4.阶跃函数的多项式逼近。

因此，给定一组适当选择的参数（α1，...，αn），PCT能够对各种信号进行高精度分析，这些信号的IF轨迹可以是时间的任何连续函数。

C.数值实验测试

拟议的PCT的性能进行评估，通过测试两个数字生成的信号。第一个信号被给出为

以200 Hz的采样频率对信号进行采样。该信号的IF轨迹是时间的非线性函数，即，f（t）=10 + 2.5t + t2/3 − t3/40（Hz）。图1A和1B中所示的信号的TFD。图5（a）-（c）和6分别由传统的STFT、CT和提出的PCT产生。在CT中，啁啾率取6π，PCT中多项式核特征参数取（α1，.，α3）=（5π，2π/3，−π/20）;由此，多项式核可以精确地跟踪信号的IF轨迹的曲线。基本上，所有三个TFD都可以呈现信号的固有时频模式。然而，很明显，由PCT产生的TFD具有最佳浓度，而由STFT产生的TFD具有最差浓度。此外，可以看出，对于CT，TFD在0-10 s处的集中度优于10-15 s处的集中度;这主要是由于信号在[0，10] s时间跨度上的IF轨迹可以很好地近似为线性啁啾信号，其频率以3 Hz/s的速度随时间增加。

图5.由（12）给出的信号的TFD（窗口大小= 512）。(a)由STFT。(b)通过CT。(c)通过PCT。

这里考虑的第二信号具有由图4所示的阶梯状函数描述的IF轨迹。STFT和PCT都用于计算该信号的TFD。当使用PCT时，多项式核特征参数如（11）中所给出。在这种情况下，PCT与一组适当的参数产生的TFD在浓度方面比STFT的TFD分析结果质量好得多，这并不奇怪。

图6.具有阶梯状IF轨迹的信号的TFD（窗口大小= 512）。(a)由STFT。(b)通过PCT。

4.结论

本文通过对传统CT的扩展，提出了一种新的时频分析方法，称为PCT。PCT不仅可以产生具有优良浓度的线性啁啾信号，其IF是时间的线性函数，而且还可以产生具有优良浓度的非线性啁啾信号，其IF是时间的非线性函数。此外，基于新提出的PCT算法，提出了一种有效的信号中频估计算法，并将该算法应用于含非线性啁啾分量、受高斯噪声严重污染的信号和转子试验台振动信号的中频估计，验证了该算法的有效性。

5.MATLAB代码及仿真结果

function [Spec,f] = Polychirplet(Sig,SampFreq,Ratio,N,WinLen);
​
% Polynomial chirplet transform
%
%     Sig       : the signal to be analyzed
%    SampFreq    : sampling frequency 
%      Ratio     : Coefficients of polynomial
%        N       : the number of frequency bins (default : length(Sig)).
%  WinLen       : the length of window used to locate the signal in time.
%
​
​
if(nargin < 3),
    error('At least 3 inputs are required!');
end
​
SigLen = length(Sig);
​
if (nargin < 5),
    WinLen = SigLen / 4;
end
​
if (nargin < 4),
    N = SigLen;
end
​
if(N > 512),
    N = 512;
end
​
RatioNum = length(Ratio);
​
dt = (0:(SigLen-1))';
dt = dt / SampFreq;
​
df = zeros(size(dt));
​
for k = 1:RatioNum,
    df = df + Ratio(k) * dt.^k;
end
​
kernel = zeros(size(dt));
​
for k = 1:RatioNum,
    kernel = kernel + Ratio(k)/(k + 1) * dt.^(k + 1);
end
​
WinLen = ceil(WinLen / 2) * 2;
t = linspace(-1,1,WinLen)';
WinFun = exp(log(0.005) * t.^2 );
WinFun = WinFun / norm(WinFun);
​
Lh = (WinLen - 1)/2; 
​
Spec = zeros(N,SigLen) ;  
Rt = zeros(N,1);
Rdt = zeros(N,1);
​
wait = waitbar(0,'Please wait...');
for iLoop = 1:SigLen,
​
    waitbar(iLoop/SigLen,wait);
    
    tau = -min([round(N/2)-1,Lh,iLoop-1]):min([round(N/2)-1,Lh,SigLen-iLoop]);
    temp = floor(iLoop + tau);
    Rt(1:length(temp)) = kernel(temp); 
    Rdt(1:length(temp)) = dt(temp);  
    
    temp1 = floor(Lh+1+tau);
    rSig = Sig(temp);
    rSig = Hilbert(real(rSig));
​
    rSig = rSig .* conj(WinFun(temp1));
    Spec(1:length(rSig),iLoop) = rSig;
    
    Spec(:,iLoop) = Spec(:,iLoop) .* exp(-j * 2.0 * pi * (Rt - df(iLoop) * Rdt)); 
​
%    ft = abs(fft(Spec(:,iLoop)));
%    Spec(:,iLoop) = ft;
    
end;
​
Spec = fft(Spec); 
Spec = abs(Spec);
​
close(wait);
​
Spec = Spec(1:round(end/2),:);
[nLevel, SigLen] = size(Spec);
​
f = [0:nLevel-1]/nLevel * SampFreq/2;
t = (0: SigLen-1)/SampFreq;
​
[fmax fmin] = FreqRange(Sig);
fmax = fmax * SampFreq;
fmin = fmin * SampFreq;
​
clf
%=====================================================%
% Plot the result                                     %
%=====================================================%
set(gcf,'Position',[20 100 350 300]);      
set(gcf,'Color','w');              
​
mesh(t,f,Spec);  
axis([min(t) max(t) fmin fmax]);
ylabel('Freq / Hz');
xlabel('Time / Sec')
Info = 'C = ';
​
for i = 1:RatioNum,
    Info = [Info,num2str(Ratio(i),4), '  '];
end
​
if RatioNum == 0,
    Info = 'C = 0';
end
​
%title(['Nonlinear Chirplet[',Info,']']);
​