【时频分析】用于TFD性能评估和比较的最优分段样条Wigner-Ville分布的设计【附MATLAB代码】

摘要

通过设计一种参考最优时频分布以及新的精度和分辨率度量，提出了一种新的时频分布性能评价方法。其动机来自于需要一种客观的、能够量化TFD的精度和分辨率、能够确定不同TFD之间的性能差异的TFD性能评估方法，并且适用于具有任意分量数目、瞬时频率和幅度的信号。通过使用分段线性调频(LFM)基对标准非平稳信号模型进行分解，并利用LFM信号的Wigner-Ville分布的最优性，提出了所提出的最优TFD，即分段Spline Wigner-Ville分布(PWWVD)。我们将所设计的PW-WVD与传统的最优TFD进行了比较，表明前者更适合作为TFD性能评估的参考。利用PW-WVD，我们推导了TFD的精度和分辨率度量，并将它们与传统方法进行了比较，分析了它们的灵敏度，形成了TFD的选择标准。我们评估了12种不同TFD的精度和分辨率，并开发了在有或没有关于信号参数的先验信息的情况下的精确TFD选择策略。结果表明，在没有关于信号参数的先验信息的情况下，紧致核分布是性能最好的TFD，并且必须根据先验信息的可用性来选择不同的TFD。

关键词：多分量信号，性能评估，分段线性调频信号，时频分布，时频测量，时频表示。

时频表示（TFR）由时频分布（TFD）获得，描述了信号随时间的频谱内容[1]。它是分析和处理谱信息随时间变化的非平稳信号的一种适当表示[1]。时频（TF）信号处理技术在诸如异常检测的各种应用中，当与仅时间或仅频率方法相比时，已经显示出上级结果分类[2]-[5]，定位[6]，[7]，频率估计[8]-[10]，信号建模与合成[11]，[12]，参数估计[13]-[15]。TFD的精度和分辨率决定了其性能[1]。高精度可确保最小伪影（交叉项），而高分辨率可增强分辨密集分量的能力[1]。目视评估和TF测量是评估TFD性能的可行工具[16]。然而，视觉评估是主观的，TF测量不能定义不同TFD之间的性能差异。因此，我们提出了一个新的TFD性能评估过程，通过设计一个参考最优的TFD和新的精度和分辨率的措施。

本文的目标是：（1）针对具有非线性瞬时频率（IF）和瞬时幅度（IA）的多分量非平稳信号，设计一种最优的时域有限差分（TFD）;（2）以所提出的分布为参考，比较和评估时域有限差分的性能。

本文的主要贡献如下：（1）设计了一个最优的分段样条Wigner-Ville分布（PW-WVD）;（2）使用所提出的分布导出TFD的精度和分辨率度量;（3）使用所提出的度量评估各种TFD的性能;（4）提供TFD选择策略，无论是否有信号参数的先验信息。

Wigner-Ville分布（WVD）产生单分量线性调频（LFM）信号的最佳TFR [1]。然而，已知在分析多分量或非线性调频信号时会遇到交叉项[17]。高性能的TFD是通过使用TF内核过滤WVD来实现的，该TF内核实现了准确度和分辨率之间的最佳权衡[1]。各种研究采用TFD来改善不同领域的传统信号处理方法的性能，例如生物医学[18]，[19]，雷达[20]，[21]和电信[22]，[23]。目前，TFD性能使用TF度量进行优化和评估，例如范数比[24]，能量集中[25]，熵[26]和特定应用标准[27]。能量集中在[28]中用于基于短时分数傅立叶变换优化频谱图。调整窗口参数和分数角以使浓度最大化。此外，在[25]中使用该度量的标准化版本来自动优化各种TFD的性能。R 'enyi熵也用于通过最小化TFD的复杂性来优化TFD的性能[29]。在[24]中，各种R ′ enyi熵和能量集中度量用于分析和量化性能。以同样的方式，特定于应用的标准，如分类精度和估计误差，通常用于评估TFD的性能[30]，[31]。

TF测量是与TFD的性能相关的代理。比较不同TFD的性能是有用的，但它不能定义这些TFD之间的性能差异。归一化瞬时分辨率（NIR）旨在减轻TF测量陷阱[32]，[33]。然而，它需要准确检测信号分量，并且在任何时刻最多两个分量具有近似相等的幅度[34]。Reinhold测度通过放松等幅限制来扩展NIR [34]。它采用自适应检测算法来识别信号分量，而不管它们的IF复杂度和幅度差。尽管如此，它仍然限于在任何时刻具有最多两个分量的信号。因此，NIR和Reinhold测量不适合分析具有任意数量分量的信号。作为比较参考的最佳分布的效用，允许直接评估的TFD的准确性和分辨率，而不依赖于代理措施。此外，参考分布适合于确定不同的TFD之间的性能差异，即使在分析不同的信号。目前，参考TFD使用理想分布[28]，[35]，基于WVD的方法[36]，参数化TFD [37]或通过匹配追踪分布（MPD）[38]，[39]导出。理想的TFD由沿着信号IF定律的刀口组成。尽管如此，它们不适合作为参考，因为它们的TF平铺方法导致不连续性。此外，基于WVD的参考分布是通过对每个信号分量的WVD求和来创建的;因此，由于交叉项的上升，它们的效用仅限于多分量LFM信号。另一方面，参数化的TFD对于特定类别的信号是最佳的。因此，它们不能推广到具有任意数量的分量和IF/IA复杂度的信号。相比之下，MPD理论上可以作为参考分布，给出具有足够多样性的字典[40]。然而，其能力受到计算复杂性和输入信号可变性的阻碍[41]。因此，当分析具有任意数量的分量和IF/IA复杂度的信号时，MPD不适合用作参考分布[42]。在单分量LFM信号的WVD最优性[1]，线性B样条理论[43]，[44]和[45]中的分段信号模型的驱动下，我们设计了具有非线性IF/IA律的多分量非平稳信号的最优TFD。此外，我们利用建议的分布作为参考，比较和评估各种TFD的性能。我们得到新的TFD的准确性和分辨率的措施，并制定精确的TFD选择策略。

本文的其余部分组织如下：第二节制定拟议PW-WVD和提出的方法来评估TFD的准确性和分辨率。之后，我们在第三节中介绍了TFD性能评估结果和开发的TFD选择策略。最后，第四节结束了工作。

Fig. 1.使用最佳PW-WVD的拟议TFD性能评估过程的流程图。工艺步骤总结如下：（1）给定一组参数，生成多分量非平稳信号z（t），（2）计算最优PW-WVD Lz（t，f），（3）计算WVD Wz（t，f），（4）优化核γ（t，f），（5）在准确度和分辨率方面评估TFD性能，（6）针对不同TFD重复步骤4-5，（7）对于不同的信号参数重复整个过程，以及（8）基于性能评估开发TFD选择策略。

所提出的使用最佳PW-WVD的TFD性能评估过程如图1所示，并总结如下：

1)给定一组参数（分量数量、IF、IA和时间支持），使用标准非平稳信号模型生成多分量信号z（t）。

2)使用分段LFM（PW-LFM）分解z（t）并计算其最优PW-WVD Lz（t，f）。

3)计算信号WVD Wz（t，f）。

4)给定一个核γ（t，f），通过最小化Lz（t，f）和TFD之间的差来优化其参数。

5)使用优化的TFD核γopt（t，f）、Wz（t，f）和Lz（t，f），使用所提出的精度α和分辨率α来评估TFD性能。

6)对不同的TFD内核重复步骤4-5，并汇总性能评估。

7)对不同的信号参数重复步骤1-6，并汇总所有性能评估。

8)当信号参数的先验信息可用或不可用时，制定TFD选择策略。

其余小节详细介绍了拟议的TFD绩效评估过程中的每个步骤。

A.信号模型

多分量非平稳信号由几个单分量组成，并形成一般情况。一个多分量有限非平稳信号z（t）被定义为：

其中，P是z（t）中的分量的数量，ap（t）和fp（t）分别是第p个分量的IA和IF，j = -1，并且p（t）是定义从tp，i到tp，f的第p个分量时间支持的矩形窗口。

B.分段线性调频分解

一个连续函数h ∈ C [0，1]可以近似，给定{tk} tk ∈ {tk}，而不需要明确知道它的真实参数形式[46]。可以通过用直线连接相邻点（tk，h（tk））和（tk+1，h（tk+1））来表示h的线性B样条逼近。当{tk}是0和1之间的无限基数点集时，这种近似成为h的分解[46]。

信号第p个分量的IF，如等式2所述。(1)，可以使用线性B样条分解如下[44]：

其中Bp，q（t）是fp（t）的第q次线性B样条。

使用等式(3)，等式中的信号定义(1)变成：

当量(5)可以用以下替代形式表示：

通过比较等式(1)等式(6)，我们观察到，求和项是等式中第p个复指数的PW-LFM分解。这种分解使我们能够利用LFM信号的WVD最优性来设计所提出的最优PW-WVD。

C.分段样条Wigner-Ville分布

给定解析信号z（t），WVD定义为[1]：

其中Kz（t，τ）是z（t）的瞬时自相关函数（IAF），Wz（t，f）是WVD，z ω（t）是z（t）的复共轭，F τ→f表示从滞后τ到频率f的傅立叶变换。众所周知，WVD遭受交叉项1。然而，它对于单分量LFM信号是最佳的。这种最优性允许我们使用信号PW-LFM分解模型来制定PW-WVD。

1)构想：使用Eq.(6)，z（t）自动项的IAF表示为：

使用公式(8)以及（10），PW-WVD被定义为：

PW-WVD是最佳的，给出了信号参数α f（t），α（t），β f（t），β f（t）。此外，它提供了一种将WVD交叉项Xz（t，f）与其自项分离的方法，即：

图二.由等式（1）中的信号模型描述的示例2分量非平稳信号的自项和交叉项。（一）.信号IA、IF和时间支持函数分别在（a-c）中演示。PW-WVD显示了（d）中的信号自动项。信号WVD及其分离的交叉项分别在（e）和（f）中示出。信号第一分量的特征在于指数衰减的IA，其形式为：exp（−0.008t），通过TF点{（20，0.4），（80，0.2），（130，0.35），（180，0.15），（230，0.3）}的四阶多项式IF，以及20秒和230秒之间的时间支持。信号第二分量的特征在于正弦IA的形式：0.4+0.3 cos（0.01πt），通过{（100，0.05），（220，0.4）}的线性IF，以及100和220秒之间的时间支持。模拟参数为：N = 256个样本，fs = 1 Hz，M = 256个样本。

图2示出了示例2分量非平稳信号的PW-WVD、WVD和WVD分离的交叉项。此外，附录B讨论了PW-WVD最优性及其性质，而附录C给出了离散PWWVD公式。

MATLAB代码

% Copyright (c) 2021 Mohammad Fathi Al-Sa'd
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the "Software"),
% to deal in the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the
% Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%
% Email: mohammad.al-sad@tuni.fi, alsad.mohamed@gmail.com
%
% The following reference should be cited whenever this script is used:
% M. Al-Sa'd, B. Boashash, M. Gabbouj, "Design of an Optimal Piece-Wise
% Spline Wigner-Ville Distribution for TFD Performance Evaluation and
% Comparison", IEEE Transactions on Signal Processing, 2021.
%
% Last Modification: 21-May-2021
%
% Description:
% This demo script produces the results that are depicted in Fig. 3 of the
% paper. It generates an example 2-component non-stationary signal and
% computes the signal ideal TFD, PW-WVD, and the MPD using chirplet atoms.
% Note that the MPD is precomputed and saved in MPD_example.mat.
​
%% Initialization
clear; close all; clc;
addpath(genpath('Functions'));
​
%% Parameters
N  = 256;
fs = 1;
t  = 0:1/fs:(N-1)/fs;
f  = 0:fs/(2*N-1):fs/2;
C_max = 128;
​
%% Generate component 1
fc = 0.3;
fd = 0.1;
fm = 0.015;
IF(1,:) = fc - fd.*sin(2*pi*fm*t);
IP(1,:) = fc*t + (fd./(2*pi*fm)).*cos(2*pi*fm*t);
W(1,:)  = ones(1,N);
IA(1,:) = ones(1,N);
z(1,:)  = IA(1,:).*exp(1j*2*pi*IP(1,:)).*W(1,:);
​
%% Generate component 2
ff = 0.025;
m  = 0.1/255;
IF(2,:) = ff + m*t;
IP(2,:) = ff*t + 0.5.*m.*t.^2;
W(2,:)  = ones(1,N);
IA(2,:) = exp(-0.0037*t);
z(2,:)  = IA(2,:).*exp(1j*2*pi*IP(2,:)).*W(2,:);
z       = sum(z);
​
%% Ideal TFD
Iz = ideal_tfd(IF,IA,W);
​
%% PW-WVD
Lz = pw_wvd(IF,IP,IA,W,fs);
​
%% MPD using chirplet atoms.
% Uncomment the 5 lines below if you want to re-compute the MPD.
% dict = chirplet_dictionary(N,fs);
% [atoms,c,idx] = mp_atoms(z,dict,1);
% Mz = mpd(atoms);
% save(['Data' filesep 'MP' filesep 'MPD_example.mat'],'Mz');
% save(['Data' filesep 'MP' filesep 'MP_example.mat'],'c','idx');
load(['Data' filesep 'MP' filesep 'MPD_example.mat']);
​
%% Ideal TFD TF-Tiling
ideal_t_tiling = repmat(t,N,1);
​
%% PW-WVD TF-Tiling
P = size(IF,1);
pw_wvd_tf_tiling = cell(1,P);
for p = 1:P
    IF_temp = IF(p,:);
    IFw = IF_temp - max(IF_temp);
    fd = 0:fs/(2*N-1):(abs(min(IFw))+fs/2);
    for i = 1:length(fd)
        pw_wvd_tf_tiling{p}(i,:) = IFw + fd(i);
    end
end
​
%% MPD TF-Tiling
load(['Data' filesep 'MP' filesep 'MP_example.mat']);
parameters = chirplet_parameters(N,fs);
mp_param = parameters(:,idx);
index = [1 4 8];
mpd_tf_tiling = cell(1,P);
for p = 1:length(index)
    IF_temp = mp_param(end,index(p)).*(t - mp_param(2,index(p)))./pi;
    IFw = IF_temp - max(IF_temp);
    fd = 0:fs/(2*N-1):(abs(min(IFw))+fs/2);
    for i = 1:length(fd)
        mpd_tf_tiling{p}(i,:) = IFw + fd(i);
    end
end
​
%% Plotting ideal TFD, PW-WVD, and MPD
x1 = 0.035; y1 = 40.5;  w1 = 0.022; h1 = 30;
x2 = 0.351;  y2 = 105.5; w2 = 0.057;  h2 = 22;
​
figure('Color',[1,1,1],'Position',[100 100 2200 550]);
f1 = subplot(1,3,1);
f1.Position = [0.125 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Iz)); axis xy; hold on;
plot(f,0.5+ideal_t_tiling(:,4:5:end),'-','linewidth',3,'Color',[0 0.4470 0.7410]);
rectangle('Position',[x1 y1 w1 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
axis([x1 x1+w1 y1 y1+h1]); caxis([0 C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.338181818181818 0.425909090909091],...
    [0.923636363636363 0.367272727272727],'Color',[0 0.4470 0.7410],'LineWidth',5);
annotation('line',[0.338181818181818 0.425909090909091],[0.16 0.28],...
    'Color',[0 0.4470 0.7410],'LineWidth',5);
f2 = subplot(1,3,2);
colormap(1-gray); imagesc(f,t,abs(Iz)); axis xy; hold on;
rectangle('Position',[x1 y1 w1 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
xlabel('Frequency (Hz)'); ylabel('Time (s)');
axis([f(1) f(end) t(1) t(end)]); caxis([0 C_max]);
set(gca,'fontweight','bold','fontsize',24,'FontName','Times');
f3 = subplot(1,3,3);
f3.Position = [0.7 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Iz)); axis xy; hold on;
plot(f,0.5+ideal_t_tiling,'-','linewidth',3,'Color',[0.4660 0.6740 0.1880]);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
axis([x2 x2+w2 y2 y2+h2]); caxis([0 C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.585454545454545 0.7],...
    [0.543636363636364 0.925454545454545],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
annotation('line',[0.585 0.700454545454545],...
    [0.476363636363636 0.158181818181818],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
set(gcf,'Units','inches'); screenposition = get(gcf,'Position');
set(gcf,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',screenposition(3:4));
​
figure('Color',[1,1,1],'Position',[100 100 2200 550]);
f1 = subplot(1,3,1);
f1.Position = [0.125 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Lz)); axis xy; hold on;
plot(0.5.*fs/(2*N-1)+pw_wvd_tf_tiling{2},t,'-','linewidth',3,'Color',[0 0.4470 0.7410]);
rectangle('Position',[x1 y1 w1 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
axis([x1 x1+w1 y1 y1+h1]); caxis([0 2*C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.338181818181818 0.425909090909091],...
    [0.923636363636363 0.367272727272727],'Color',[0 0.4470 0.7410],'LineWidth',5);
annotation('line',[0.338181818181818 0.425909090909091],[0.16 0.28],...
    'Color',[0 0.4470 0.7410],'LineWidth',5);
f2 = subplot(1,3,2);
colormap(1-gray); imagesc(f,t,abs(Lz)); axis xy; hold on;
rectangle('Position',[x1 y1 w1 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
xlabel('Frequency (Hz)'); ylabel('Time (s)');
axis([f(1) f(end) t(1) t(end)]); caxis([0 2*C_max]);
set(gca,'fontweight','bold','fontsize',24,'FontName','Times');
f3 = subplot(1,3,3);
f3.Position = [0.7 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Lz)); axis xy; hold on;
plot(0.5.*fs/(2*N-1)+pw_wvd_tf_tiling{1}(1:2:end,:),t,'-','linewidth',3,'Color',[0.4660 0.6740 0.1880]);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
axis([x2 x2+w2 y2 y2+h2]); caxis([0 2*C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.585454545454545 0.7],...
    [0.543636363636364 0.925454545454545],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
annotation('line',[0.585 0.700454545454545],...
    [0.476363636363636 0.158181818181818],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
set(gcf,'Units','inches'); screenposition = get(gcf,'Position');
set(gcf,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',screenposition(3:4));
​
figure('Color',[1,1,1],'Position',[100 100 2200 550]);
f1 = subplot(1,3,1);
f1.Position = [0.125 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Mz)); axis xy; hold on;
plot(0.5.*fs/(2*N-1)+mpd_tf_tiling{1},t,'-','linewidth',3,'Color',[0 0.4470 0.7410]);
rectangle('Position',[x1+0.002 y1 w1-0.002 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
axis([x1+0.002 x1+w1 y1 y1+h1]); caxis([0 C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.338181818181818 0.425909090909091],...
    [0.923636363636363 0.367272727272727],'Color',[0 0.4470 0.7410],'LineWidth',5);
annotation('line',[0.338181818181818 0.425909090909091],[0.16 0.28],...
    'Color',[0 0.4470 0.7410],'LineWidth',5);
f2 = subplot(1,3,2);
colormap(1-gray); imagesc(f,t,abs(Mz)); axis xy; hold on;
rectangle('Position',[x1 y1 w1 h1],'edgecolor',[0 0.4470 0.7410],'linewidth',5);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
xlabel('Frequency (Hz)'); ylabel('Time (s)');
axis([f(1) f(end) t(1) t(end)]); caxis([0 C_max]);
set(gca,'fontweight','bold','fontsize',24,'FontName','Times');
f3 = subplot(1,3,3);
f3.Position = [0.7 0.1589 0.2134 0.7661];
colormap(1-gray); imagesc(f,t,abs(Mz)); axis xy; hold on;
plot(0.5.*fs/(2*N-1)+mpd_tf_tiling{2}(1:5:end,1:117),t(1:117),'-','linewidth',3,'Color',[0.4660 0.6740 0.1880]);
plot(0.004+mpd_tf_tiling{3}(1:5:end,117:end),t(117:end),'-','linewidth',3,'Color',[0.4660 0.6740 0.1880]);
rectangle('Position',[x2 y2 w2 h2],'edgecolor',[0.4660 0.6740 0.1880],'linewidth',5);
axis([x2 x2+w2 y2 y2+h2]); caxis([0 0.3.*C_max]);
set(gca,'XtickLabel','','YtickLabel','');
annotation('line',[0.585454545454545 0.7],...
    [0.543636363636364 0.925454545454545],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
annotation('line',[0.585 0.700454545454545],...
    [0.476363636363636 0.158181818181818],'Color',[0.4660 0.6740 0.1880],'LineWidth',5);
set(gcf,'Units','inches'); screenposition = get(gcf,'Position');
set(gcf,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',screenposition(3:4));
​
%% Saving
opt = input('Do you want tosaveall results (Y/N)\n','s');
if(opt == 'y' || opt == 'Y')
    print(1,'comparison_Iz','-dpdf','-r400');
    print(2,'comparison_Lz','-dpdf','-r400');
    print(3,'comparison_Mz','-dpdf','-r400');
end