数字信号处理-时域离散随机信号处理课件：时频分析

时频分析5.1引言

5.2短时傅里叶变换

5.3维格纳变换(WD)5.4时域离散信号的维格纳变换

5.5时频分布的统一表示式

5.6时频分析在编队目标架次检测中的应用

5.1引

言

传统的信号分析与处理的数学工具是傅里叶变换，它的正变换和逆变换分别用下面两式表示：

(5.1.1)(5.1.2)式中ω是一个连续变量,限制了用计算机在频域进行分析与处理，而离散傅里叶变换(DFT)将频域离散化,使之借助计算机可以在时域也可以在频域对信号进行分析与处理。由于傅里叶变换物理概念清晰,同时也是正交变换,因此长期以来科技界及各工程领域广泛使用傅里叶变换和离散傅里叶变换。

X(ejω)称为信号x(n)的频谱，它表示了信号在频域的分布规律。也可以用下面公式表示:（5.1.3）

e(ω)称为信号x(n)的能量谱，它仅包含信号的幅度信息。但对于能量无限信号，如周期信号、平稳随机信号等，傅里叶变换不存在，可以用功率谱密度（简称功率谱）P(ejω)表示:（5.1.4）

式中rxx(m)是x(n)的自相关函数。频谱、能量谱以及功率谱都是信号变换到频域的一种表示方法，对于频谱不随时间变化的确定性信号以及平稳随机信号都可以用它们进行分析和处理。5.2短时傅里叶变换

5.2.1短时傅里叶变换的定义及其物理解释

1.短时傅里叶变换的定义短时傅里叶变换的定义有两种形式,下面分别叙述。

(1)定义一:(5.2.1)

式中w(n)是一个窗函数，其作用是取出x(n)在n时刻附近的一小段信号进行傅里叶变换，当n变化时，窗函数随n移动，从而得到信号频谱随时间n变化的规律。此时的傅里叶变换是一个二维域(n,ω)的函数。窗函数沿时间轴移动情况如图5.2.1所示。图

5.2.1窗函数的移动

令n′=n-m

将n′代入定义一中,再将n′用m代替,可得到第二种定义形式。

(2)定义二：

2.短时傅里叶变换的物理解释

对以上STFT的定义形式，从傅里叶变换和线性滤波两个角度，可以有两种不同的物理解释。

(1)由傅里叶变换角度解释。按照(5.2.1)式，STFT可以看作n是参变量，x(m)w(n-m)对m的傅里叶变换，因此它是（n,ω）的函数。因为STFT是x(m)w(n-m)的傅里叶变换，可以用x(m)和w(n-m)分别的傅里叶变换的卷积表示。设：那么

如果再将θ改换成-θ，

得到

(5.2.3)上式是STFT定义的一种频域表示形式。这里如果x(n)是时变信号，式中用了它的傅里叶变换，是不合适的，但可以理解为信号在时间窗外变为0以后，取信号的傅里叶变换;或者说是时间窗内的信号傅里叶变换的平滑形式。（2）

由线性滤波角度解释。将定义一重写如下：

上式表明,短时傅里叶变换可以看成x(n)e-jωn与w(n)的线性卷积，如将w(n)看成一个低通滤波器的单位脉冲响应，短时傅里叶变换则可用图5.2.2表示。图5.2.2表明,首先将信号x(n)调制到-ω，然后通过低通滤波器w(n)，其输出就是短时傅里叶变换。实质上是将x(n)在ω附近的频谱搬移到零频处，作为短时傅里叶变换。为使其频率分辨率高，希望w(n)是一个低通窄带滤波器，带外衰减愈大愈好。

利用定义二可以得到线性滤波的另一种物理解释，

将定义二重写如下:公式中求和号部分可看成w(n)ejωn与x(n)的线性卷积，因此上式可以写成式中w(n)是低通滤波器，w(n)ejωn就是以ω为中心的带通滤波器。按照上式，STFT就是信号首先通过带通滤波器，选出以ω为中心的频谱，再乘以exp(-jωn)，将选出的频谱搬移到零频处。

短时傅里叶变换如按照定义二的物理解释，则可用图

5.2.3表示。

图5.2.2定义一的物理解释图

5.2.3定义二的物理解释

5.2.2短时傅里叶变换的性质短时傅里叶变换是建立在一般傅里叶变换基础上的一种变换，因此它具有许多和傅里叶变换相似的性质。

1.线性性质

设

z(n)=c·x(n)+d·y(n),c,d

为常数,则

(5.2.4)2.频移性质——调制特性设 ,则

(5.2.5)3.时移特性设x(n)=y(n-n0),则

(5.2.6)证明

以上说明STFT具有频移不变性，但不具有时移不变性，相差一个相位因子。

4.共轭对称性

当信号是实信号时,短时傅里叶变换和一般傅里叶变换一样具有共轭对称性,即

(5.2.7)因此,其实部是偶函数，虚部是奇函数。

5.由短时傅里叶变换恢复信号由定义(5.2.1)式得到短时傅里叶变换的反变换为

设n=m,则

(5.2.8)

只要w(0)≠0,可以由STFTX(n,ω)准确地恢复信号x(n)。

5.2.3短时傅里叶变换的时间、频率分辨率由定义可知，STFT实际分析的是信号的局部谱，局部谱的特性决定于该局部内的信号，也决定于窗函数的形状和长度。为了了解窗函数的影响，假设窗函数取两种极端情况。第一种极端情况是取w(n)=1，

-∞<n<∞,此时

这种情况下，STFT退化为信号的傅里叶变换，没有任何时间分辨率，却有最好的频域分辨率。第二种极端情况是取w(n)=δ(n)，

此时

STFT退化为信号，有理想的时间分辨率，但不提供任何频率分辨率。

短时傅里叶变换由于使用了一个可移动的时间窗函数，使其具有一定的时间分辨率。显然,短时傅里叶变换的时间分辨率取决于窗函数w(n)的长度。为了提高信号的时间分辨率,希望w(n)的长度愈短愈好。但频域分辨率取决于w(n)窗函数的频域函数宽度，也就是低通滤波器w(n)的带宽或者说带通滤波器w(n)ejωn的带宽，为了提高频域分辨率，希望尽量加宽w(n)窗口宽度，这样必然又会降低时域分辨率。因此,

STFT的时间分辨率和频率分辨率不能同时任意提高。这种时域分辨率和频域分辨率相互制约的性质，也正反映了已为理论所证明了的“不确定原理”：（5.2.9）

式中Δt表示信号有效持续时间，Δf表示信号的有效带宽。上面公式说明,对于窗函数，它的时间宽度和在频率域的宽度不能同时任意小,也就是说,频域分辨率和时域分辨率不能同时任意小。但可以选择合适的窗函数，使Δt和Δf都比较小,其乘积接近于1/(4π)。窗函数的形式有很多,可以证明从有效时宽和有效频宽乘积为最小的意义上讲，高斯波形信号是最好的，但是它在时间轴和频率轴上是无限扩张的，因此它并不是一种最好的波形。我们知道，不可能存在既是带限又是时限的信号波形，实际应用中采用放松条件，研究在有限时宽的情况下，使频率有效带宽为最小的波形是什么，或者研究在有限带宽情况下，使时宽最小的波形是什么，这部分内容可参考文献［8］、［2］。5.2.4短时傅里叶变换的实现

1.短时傅里叶变换的时域、频域采样

我们已经知道,短时傅里叶变换是低通滤波器w(n)的输出，假设其有效带宽为Δ，对应的模拟滤波器的有效带宽为B,式中fs是x(n)的采样频率，那么对应该模拟滤波器的时域采样频率应是带宽的两倍(2B)以上,最小采样频率为2B即是最小再次采样率。上式表明该采样率是信号采样率fs的P倍,也就是说,二次采样间隔最大为1/P，即窗口每次移动的最大间隔是1/P。例如:窗函数选用哈明窗,长度为L,带宽近似为4π/L，设:fs=10kHz,L=100,则

即最大采样间隔为25，n的取值为0,25,50,75,…。以上计算的再次采样率是最小采样率，采样间隔是最大采样间隔。对于频率域采样，假设在周期2π中采样M点，为不发生时域混叠，

要求M≥L。

2.用FFT计算STFT假设在频率域等间隔采样M点，

k=0,1,2,3,…,M-1k=0,1,2,3,…,M-1按照定义一，有

（5.2.10）

令m=l+n，

上式变为

（5.2.11）

上式的求和区间是(-∞,∞)，可以按照长度为M的区间进行划分。一个个区间计算后，

再求和，这样上式变成为

（5.2.12）

令l=m+rM,同时考虑到 ，得到(5.2.13)

式中

(5.2.14)m=0,1,2,3,…,M-1在(5.2.13)式中，对任何固定n

值，求和项可以用M点的FFT进行计算，其中信号用（5.2.14）式计算。根据（5.2.13）式和(5.2.14)式,由x(n)计算STFT的过程如图5.2.4所示。图

5.2.4用FFT计算短时傅里叶变换3.用滤波器组法实现短时傅里叶变换假设在频率域采样M点,采样点的频率为

k=0,1,2,3,…,M-1将ωk代入定义(5.2.2)式中，得到

令

(5.2.15)则

(5.2.16)令

(5.2.17)则

(5.2.18)这样对应M个采样点频率,形成M个通道。(5.2.15)式即是每个通道的带通滤波器的单位取样响应,(5.2.17)式即是每个带通滤波器的输出,（5.2.18）式表示每个通道的STFT输出。

图

5.2.5

STFT一个通道的原理框图

4.短时傅里叶变换的综合

由短时傅里叶变换恢复时域信号称为综合。下面先推导短时傅里叶变换M个通道的等效传输函数。将M个带通滤波器的输出相加，输出用y(n)表示,从x(n)到y(n)等效单位取样响应用h(n)表示,(5.2.19)(5.2.20)式中

等效传输函数用H(ejω)表示，即是上式的傅里叶变换。

(5.12.21)式中

(5.12.22)y(n)=x(n)*h(n)(5.2.23)(5.2.24)(5.2.23)式表示等效传输函数是M个带通滤波器传输函数的叠加。下面证明公式：

M≥L

(5.2.25)证明

对于一定的ω，令ω=ωθ,ωk’=ωθ-ωk带入上式左边，得到上式最右边求和号是频域信号 的IDFT，即w(-n)。

另外,我们知道频域采样，使时域以采样点数M为周期进行周期性延拓，

因此

n=0,1,2,3,…,M-1这里M≥L

保证不出现时域混叠,令n=0,r=0,带入上式，最后得到

(证毕)将上式代入(5.2.23)式,得到

H(ejω)=Mw(0),h(n)=F-1T［H(ejω)］=Mw(0)δ(n)

y(n)=x(n)*h(n)=Mw(0)x(n)(5.2.26)(5.2.26)式说明当M≥L时,用短时傅里叶变换可以恢复原时域信号,恢复的可能性与窗函数的形状无关。但是窗函数的形状会影响其时间与频率分辨率。由(5.2.18)式得到

(5.2.27)按照(5.2.19)式，

有

(5.2.28)上式表明，由STFT恢复原信号，也就是将每一路的STFT输出乘以exp(jωkn)，然后进行相加得到时域信号y(n)。下面将x(n)的短时傅里叶分析和y(n)信号的恢复画在一起，如图5.2.6所示。图

5.2.6STFT分析与综合的原理图

我们知道,STFT的输出的带宽决定于低通滤波器的带宽，其带宽一般比输入信号的带宽小得多，如果输入信号的采样率满足采样定理，那么在STFT的输出完全可以降低采样率，对STFT的输出进行二次采样，减少数据量，即运算量。为此,在信号STFT输出端加一个抽取器，信号综合时，再进行插值，以恢复原来的采样率。这样,采用滤波器组实现STFT分析与综合的原理框图如图5.2.7所示,图中(a)与(b)分别对应定义一和定义二。二次采样率用输入信号的带宽与低通滤波器的带宽之比进行计算，设输入信号带宽用Δx表示，低通滤波器的带宽用Δω表示，图中的L用下式计算：

(5.2.29)

上面我们研究了STFT的计算方法，这里主要的问题是计算的有效性和时间窗的选择。为了提高STFT计算的有效性，一些学者研究了递归算法，给出了矩形滑动窗递推算法和矩形扩展窗递推算法（矩形窗的宽度作为递归形式进行扩展）;为了避免矩形窗带来的截断效应，又给出哈明（Hamming）窗和汉宁（Hanning）窗的改进递归算法，这样又增加了计算的复杂性。文献［3］给出了一种新的递归算法，它采用了全极点的滑动窗，它不仅克服了采用矩形窗带来的截断效应，还比用哈明窗和汉宁窗的计算量和存储量少得多。这方面内容请参考文献［2］、［3］。

图

5.2.7用滤波器组实现STFT的一个通道

5.3维格纳变换(WD)5.3.1WD的定义

确定性时间连续信号的WD定义如下:(5.3.1)定义表明这种变换是把过去某一时间信号乘上未来某一时间信号，再对两个信号时间差τ求傅里叶变换得到。

令

(5.3.2)将rxx(t,τ)称为瞬时自相关函数，那么WD就是信号瞬时自相关函数的傅里叶变换。x(t)在频率域的WD分布定义如下:(5.3.3)对于两个连续时间信号x(t)与y(t),互WD定义为

(5.3.4)同样,它们在频率域的互WD定义如下:(5.3.5)式中,X(Ω)是x(t)的傅里叶变换;Y(Ω)是y(t)的傅里叶变换。(5.3.1)、（5.3.4）式是在时域的定义形式，(5.3.3)、（5.3.5）式是在频域的定义形式。可以证明,时域和频域的定义形式有下面关系：(5.3.6)(5.3.7)5.3.2WD的性质

1.WD的实数性和对称性

(1)WD是t和Ω的实函数。

证明

对(5.3.1）式的两边取共轭,得到

令τ′=-τ,则

因此

(5.3.8)(2)如果x(t)是实信号，则WD是频率的偶函数。

即

(5.3.9)证明

将定义中的Ω换成-Ω,得到

因此

(3)对于互WD,则具有如下性质:(5.3.10)

2.边缘积分特性

(1)在固定时刻t，WD沿全频率轴的积分等于在t时刻信号的瞬时功率Px(t)(也称时间边缘特性)，即（5.3.11）

证明

由定义(5.3.1)式,得到

令

那么

令t1=t2=t,τ=0,则得到

(2)在固定频率Ω,WD沿全时间轴的积分等于该频率的能量密度Px(Ω)(也称频率边缘特性)，即(5.3.12)证明

由定义(5.3.3)式，得到

令

,那么

代入上式，

得到

令Ω1=Ω2=Ω,因此

(3)WD分布在整个（t，Ω）平面上,对t，Ω的双重积分等于信号的总能量E，即

(5.3.13)利用(5.3.11)式和(5.3.12)式,可以推出

(5.3.14)

用以上性质可对WD进行能量化解释,时间边缘特性为信号的瞬时功率,频率边缘特性为信号的能谱密度,总能量E将Px(t)和Px(Ω)联系起来,因此WD是一种能量化的时频表示,但不能把WD解释为在时间—频率上每一点的时频能量密度,因为WD有时可能是负的,且由不确定原理不允许在某个特定的时间—频率处有能量这一概念。

3.WD的运算性质

(1)时移与频移的不变性。将x(t)时移τ，相应的其WD分布也时移τ，用公式表示如下:

如果x(t)=y(t-τ)，那么（5.3.15）

x(t)用调制，其WD分布也频移Ω0，用公式表示如下:

如果 ,那么（5.3.16）

(2)两信号的时域卷积等于两信号分别的WD在时间轴上的卷积，即:

如果y(t)=x(t)*h(t),则

（5.3.17）

(3)如果两信号的相乘，它们的傅里叶变换服从卷积关系，则和它们对应的WD在频率轴上也服从卷积关系，用公式表示如下:

如果y(t)=x(t)h(t),则

如果x(t)表示信号，h(t)表示窗函数，此性质表明信号加窗处理时，只影响频率分辨率，不影响时间分辨率。（5.3.18）

(4)两信号相加,设z(t)=x(t)+y(t),则

式中第三项称为交叉(干扰)项，

其性质意义在后面介绍。

（5.3.19）

4.WD的时限性和带限性——区域性信号的维格纳分布的时宽与频宽，与信号本身的时宽与频宽相同，即:若x(t)限制在t1≤t≤t2中，则它的WD分布也限制在同一时间域中；若x(t)的傅里叶变换X(jΩ)限制在Ω1≤Ω≤Ω2中，则它的WD分布也限制在同一频率域中，用公式表示如下：如果

则

（5.3.20）

如果

则

（5.3.21）

利用该性质，又可以得到下面结论：（1）

因果信号x(t)的WD也是因果的，即:若

t≥0t<0则

（5.3.22）

（2）解析信号z(t)的傅里叶变换限制在频率的正半轴，解析信号z(t)的WD也限制在Ω≥0的上半平面，

即

（5.3.23）

5.可逆性

由定义得到

令

则

代入上式，

得到

令t2=0,再将t1用t代替，得到

上式说明,信号可以由其WD分布进行重建，仅仅缺少初相位信息。

(5.3.24)

5.3.3常用信号的WD举例 例

5.3.1求其WD。

解

下面确定对τ的积分限：

因此

|t|<T

|t|>T

(5.3.25)

上式表明WD在时间轴上限制在-T～T之间，在频率轴上是sinx/x形式，最大值在（T,0）处。

波形图如图5.3.1所示。

图

5.3.1例

5.3.1图

例5.3.2 ,求其WD。解

该例题的信号是一个复正弦信号，可以看作平稳随机信号，其WD分布与时间无关，对任意时间都是一个在Ω=Ω0处的δ函数，

如图

5.3.2所示。

(5.3.26)图

5.3.2例

5.3.2图

例5.3.3

已知x(t)=Acos(Ω0t),求其WD。

解按照例5.3.2和(5.3.19)式，推导如下:上面结果中的第一项表示信号，第二项即是交叉干扰项。交叉干扰项将在后面介绍。x(t)的WD如图

5.3.3所示。

图

5.3.3例

5.3.3图

例5.3.4

已知 ,求其WD。解

(5.3.28)

x(t)是一个线性调频信号，其WD清楚地表示出功率谱随时间线性变化的性质。WD如图

5.3.4所示。

图

5.3.4例

5.3.4图

例5.3.5

已知 ,求互WD。

解

(5.3.29)例5.3.6已知高斯信号,求其WD。解

上式表明，

高斯信号的WD在时间上和频率上有相同的波形。

(5.3.30)5.3.4关于二次时频分布中的交叉(干扰)项前面曾说过WD是在时间、频率二维域中的能量分布函数，必然是信号的二次型，或者说是一种双线性变换。这里双线性变换指的是x(t+τ/2)x*(t-τ/2)的傅里叶变换。这样就使能量分布函数不服从线性叠加原理。我们知道,短时傅里叶变换是一种线性变换，它可以表示信号频谱随时间变化的规律。如果取它的模的平方，则可以粗略表示信号在时间、频率二维域中的能量分布。一般把短时傅里叶变换模的平方称为谱图。谱图是一种能量分布函数，也不服从线性叠加原理，两个信号之和的谱图并不等于它们分别的谱图的和，还存在第三项即交叉项。下面推导说明。假设信号x(t)谱图用SPECx(t,Ω)表示,（5.3.31）

设

式中,a,b为实常数。

式中

(5.3.32)(5.3.32)式表明,两信号相加的谱图并不等于两个信号的谱图之和，其谱图包括三部分，第一部分和第二部分分别是两个信号的谱图，第三部分即是交叉项。交叉项的最大幅度是两个信号谱图幅度乘积的两倍，而且交叉项幅度受到一个余弦函数的调制，其相位变化服从两个谱图相位差的变化规律。如果两个信号自身的谱图没有重叠部分，则第三部分为零，即没有交叉项，

这是谱图的优点。

事实上,任何二次时频分布都不服从线性叠加原理，而服从二次叠加原理。下面介绍二次叠加原理。信号x(t)的二次时频分布用Tx(t,Ω)表示，设

（5.3.33）

则

（5.3.34）

式中：和分别称为x1(t)和x2(t)的自时谱;和 分别称为x1(t)对x2(t)和x2(t)对x1(t)的互时谱。这种互时谱形成了二次时频分布的交叉项。

下面再分析WD中的交叉项，

将(5.3.33)式带入信号的双线性变换中，

得到

(5.3.35)式中，第一、二项形成自时谱，第三、四项则形成互时谱，即交叉项。对照(5.3.34)式，WD服从二次叠加原理。对于有N个分量的信号，

二次叠加原理用下式表示：

设

,则

（5.3.36）

k≠l

例5.3.7已知 ,式中c1,c2为实数，求其WD。解

(5.3.37)

该例表明，两个复单频信号频谱不重叠（Ω1≠Ω2）时，它的时频分布除了自时谱以外，仍有交叉项;交叉项处在两个信号频率连线的中点，幅度受到两信号频率差值的余弦波调制，不能保证时频分布是非负函数，其最大幅度是两信号幅度乘积的4π倍。

例5.3.3是一个高频实信号，由于分成了两个复指数信号，即存在负频率分量，在零频率处形成交叉干扰项。因解析信号没有负频率分量，先将实信号转变成解析信号，再进行WD分析，可消除这种频谱正负部分之间的交叉干扰项。下面对解析信号作简单介绍。假设x(t)是实的连续时间信号，用z(t)表示对应的解析信号。

z(t)定义为

（5.3.38）

式中是x(t)的Hilbert变换，或者说是x(t)通过一个Hilbert变换器形成的。Hilbert变换器的传输函数为(5.3.39)或者

H(jΩ)=-jsgn(Ω)(5.3.40)式中

(5.3.41)Hilbert变换器的单位冲激响应为

(5.3.42)x(t)通过Hilbert变换器后，其输出为

(5.3.43)这样解析信号z(t)和x(t)之间的关系为

(5.3.44)将上式进行傅里叶变换，

得到

Z(jΩ)=X(jΩ)+jH(jΩ)X(jΩ)=X(jΩ)［1+jH(jΩ)］

(5.3.45)将（5.3.39）式带入上式，

得到

(5.3.46)上式表明,解析信号的频谱只分布在正频率范围，是由实信号频谱的正的部分乘以2构成的;负频率部分为0。

用上式可以求解析信号。例如ejΩt，其频谱是在Ω处的δ函数，如果Ω是负的，那么也就没有正频率存在，因此它的解析信号为

Ω<0Ω>0例

5.3.8

求

f0≠0的解析信号。

解观察（5.3.39）式，它的Hilbert变换当Ω>0时，需将信号的相位变化-90°，幅度不变，因此 ,那么，x(t)的解析信号为

或者由

取其正频率部分的两倍，得到的解析信号和上式一样。

5.4时域离散信号的维格纳变换

5.4.1时域离散信号的WD定义

1.Classen提出的定义

按照连续时间信号维格纳变换定义（5.3.1）式，可以引申出Classen提出的时域离散信号的WD定义:（5.4.1）

令k′=2k，

得到

（5.4.2）

注意式中的ω是数字频率，但频域的重复周期不是2π，而是π。因此下面公式成立:

（5.4.3）

对应的频域定义为

（5.4.4）

Classen提出的定义应用较广，保持了连续时间信号WD定义中有关时间上的一些特性，但频域上的一些特性被破坏了，

例如,时间域的边缘特性为

（5.4.5）

它不同于连续时间信号WD在时间上的积分等于信号在某一频率的能量密度((5.3.12)式)。

2.Peyrin提出的定义将（5.4.1）式中的系数2去掉，

写成下式:(5.4.6)令:n+m′/2=k,n-m′/2=n′-k,则k=(n′+m′)/2,这样(5.4.6)式变成

令:n+m′/2=k,n-m′/2=n′-k,则k=(n′+m′)/2,这样(5.4.6)式变成

(5.4.7)上式就是Piyrin提出的定义。它与Classen提出的定义相比，下面关系式成立:（5.4.8）

式中,上标c代表Classen;上标p代表Peyrin。该式表明等式左面的点只是等式右边时间变量为偶数点的结果。因此,Peyrin提出的定义中包含有更多的信息，其中时间变量为奇数点的信息是Classen定义中缺少的。Peyrin提出的定义还有其它性质和优点，请参考文献［2］。

5.4.2利用FFT计算维格纳分布维格纳分布的计算量很大，目前各种快速算法还不能从根本上解决实时处理的问题。下面介绍用FFT计算维格纳分布的方法。如果用离散哈特莱变换（DHT）计算，计算的复杂性可由三倍复FFT减少到三倍实FFT的计算量，这部分内容可参考文献［2］。

离散时间信号的WD为

(5.4.9)

为用FFT进行计算，需对信号进行加窗处理，并且将频率域离散化。假设窗函数用w(l)表示，它的时宽为2L-1。且当|l|≥L时，

w(l)=0。加窗后的WD为

(5.4.10)

WDx(n,ω)的频域周期是π，若在一个周期内采样N点(N=2L-1)，采样间隔为Δω=π/N，为便于计算，再在尾部加个

0，使N=2L，令

(5.4.11)由于FFT的计算域是正的（l=0,1,…,N-1)，重新排列序列如下:l=0,1,…,L-1l=L,…,2L-1(5.4.12)最后得到

(5.4.13)上式就是用FFT计算WD的公式。

5.5时频分布的统一表示式5.5.1模糊函数及其和WD之间的关系模糊函数也是一种常用的时频表示，它广泛应用于雷达、水声等领域。本节主要介绍它的定义及其和WD之间的关系。

WD分布是对信号的双线性变换x(t+τ/2)x*(t-τ/2)关于τ作傅里叶变换，如果对该双线性变换关于时间t作傅里叶变换，则得到模糊函数的定义，公式为（5.5.1）

另外，对应WD的频率域定义（5.3.3）式，模糊函数在频率域的定义是

（5.5.2）

而且

（5.5.3）

式中,X(f)是x(t)的傅里叶变换;t为时间，τ为时延;f为频率;θ为频偏。为分析模糊函数和WD之间的关系，定义:（5.5.4）

（5.5.5）

rx(t,τ)称为瞬时自相关函数，Rx(f,θ)称为瞬间频自相关函数。这样,WD是瞬时自相关函数关于τ的傅里叶变换，模糊函数是瞬时自相关函数关于时间t的傅里叶变换。按照定义，模糊函数和WD还可以表示成以下各式：

（5.5.6）

（5.5.7）

（5.5.8）

（5.5.9）

按照定义，可以证明模糊函数具有以下性质:

（1）

时移性。令x(t)=y(t-t0),则

（5.5.10）

(2)频移性。令

,则

（5.5.11）

(3)滤波。令

（5.5.12）

则

（5.5.13）

(4)调制。令y(t)=x(t)m(t),则

(5.5.14)下面推导信号的WD和模糊函数之间的关系。

将(5.5.6)式重写如下:(5.5.15)对上式进行傅里叶反变换，得到瞬时自相关函数和WD之间的关系:(5.5.16)再对上式t作傅里叶变换，得到模糊函数和WD之间的关系:（5.5.17）

按照(5.5.8)式，还可以得到下式:（5.5.18）

(5.5.17)式和(5.5.18)式表明了WD和模糊函数之间存在如公式那样的二维傅氏变换关系。另外,还可以推导出瞬时自相关和瞬间频自相关之间的关系，如下式:(5.5.19)(5.5.20)

总结以上，(5.5.6)式～(5.5.9)式和(5.5.17)式～(5.5.20)式，将rx(t,τ)、WDx(t,f)、Rx(θ,τ)以及Ax(θ,τ)联系起来，如图5.5.1所示。图中 表示将τ映射为f的傅氏变换，反过来,则表示由f映射为τ的傅氏反变换，其它类似。图5.5.1rx(t,τ)、WDx(t,f)、Rx(θ,f)、Ax(θ,τ)之间的关系

这四个函数有四个变量，即时间变量t、时间延迟τ、频率f、频偏θ，共形成了四个域，即：（1）时频域(t,f)，对应WDx(t,f);

（2）瞬时相关域（t,τ），对应rx(t,τ);

（3）谱相关域（θ,f），对应Rx(f,θ);

（4）模糊域（θ,τ），对应Ax(τ,θ)。WDx(t,f)和Ax(τ,θ)是信号的两个不同的时频表示，按照上述分析，还应注意它们下面的两个不同点：（1）WD是能量化的时频表示，存在时间边缘特性Px(t)和频率边缘特性Px(Ω)，公式重写如下:(5.5.21)(5.5.22)Px(t)称为信号的瞬时功率，Px(f)称为信号的能谱密度，信号的总能量为

(5.5.23)模糊函数是相关化的时频表示，将模糊函数的定义重写如下:(5.5.24)(5.5.25)在(5.5.24)式中,令θ=0,得到时间相关化边缘特性rx(τ),(5.5.26)在(5.5.25)式中,令τ=0,得到频率相关化边缘特性Rx(θ),rx(τ)称为瞬时相关,Rx(θ)称为谱相关，因此模糊函数将瞬时相关和谱相关两个概念综合在一起。

（2）WD满足时频移不变性质，即满足（5.3.15）式和(5.3.16)式,或者用下式表示:如果 ,则

(5.5.28)而模糊函数满足相关化移不变性质,用公式表示如下:

如果 ,则

(5.5.29)相关化移不变性质来源于瞬时相关和谱相关移位性质得到的术语,用公式表示如下:如果 ,则

(5.5.30)总结起来,WD和模糊函数是信号的两种不同的时频表示,性质不同,相互关系是(5.5.17)、(5.5.18)式表示的二维傅氏变换对。这两种时频表示的一些对偶关系如表5.5.1所示。

表

5.5.1WD和模糊函数的一些对偶关系

5.5.2Cohen类时频分布

维格纳分布和模糊函数是实际中比较常用的两种时频分析,此外还有许多种时频分析,尽管它们的形式和性质不尽相同,但彼此常有一定的联系和共同点,因此可用一个统一的表达式表示出来,这就是L.Cohen提出的广义双线性时频表示,公式为

(5.5.31)式中φ(θ,τ)表示核函数。采用不同的核函数可以得到不同的时频分布，时频分布的各种性质要求,则反映在对核函数的约束条件上。把满足（5.5.31）式的时频分布，统称为Cohen类时频分布。

WD时频分布是Cohen类时频分布中最基本的时频分布。下面证明当核函数φ(θ,τ)=1时，

Cohen类时频分布将转换成WD。

将φ(θ,τ)=1带入（5.5.31）式,得到

下面介绍时频分布的性质对核函数要求的约束条件：（1）如要求Cx(t,Ω)=C*x(t,Ω)，即要求时频分布是实的，充要条件是要求核函数满足下式:

（2）

如要求时频分布具有时间边缘特性和频率边缘特性，

即若要求：

则要求：

φ(θ,0)=1若要求

则要求：

（3）

如要求:

则要求：

（4）

如要求时移不变和频移不变，

即若要求：

则要求φ(θ,τ)独立于t和f。

（5）

如能抑制Cx(t,f)的交叉项，则要求φ(θ,τ)为低通滤波。

（6）如Cx(t,f)具有时间支持特性，即当|t|>ta时，x（t）=0，Cx(t,f)=0,则要求:

（7）如Cx(t,f)具有频率支持特性，即当|f|>fc时,X(f)=0,Cx(t,f)=0,则要求:|θ|<2|fc|

下面介绍Cohen类时频分布的四种表示形式。

(1)二维滤波的维格纳分布。如果对核函数φ(θ,τ)作二维傅里叶变换,得到

(5.5.32)Φ(t,f)称为时频域核函数,φ(θ,τ)称为模糊域核函数。由上式得到

(5.5.33)又由WD的定义得到

(5.5.34)将(5.5.33)式和(5.5.34)式带入(5.5.31)式,得到

(5.5.35)上式就是维格纳分布与时频域核函数Φ(t,f)的二维褶积,也称为广义维格纳分布,因此Cx(t,f)类时频分布可以理解为二维滤波的维格纳分布。(2)广义模糊函数M(θ,τ)的M(-θ,τ)的二维傅里叶变换。

(5.5.36)

式中

(5.5.37)将上式带入(5.5.31)式,得到(5.5.36)式。M（θ,τ）称为广义模糊函数,它是由模糊域核函数φ(θ,τ)对Ax(θ,τ)加权所得的模糊函数。因此,(5.5.36)式可理解为广义模糊函数M(θ，τ)的M(-θ,τ)的二维傅里叶变换。

(3)广义自相关函数对时延τ的一维傅里叶变换。

(5.5.38)式中

(5.5.39)rx′(t,τ)称为广义自相关函数,而

(5.5.40)ψ(t,τ)称为瞬时相关域的核函数。

由(5.4.40)式得到

(5.5.41)(5.5.39)式和(5.5.41)式表明广义自相关函数为信号的瞬时自相关函数rx(t,τ)与瞬时相关域的核函数ψ(t,τ)对t的褶积,因此(5.5.38)式可以理解为广义自相关函数rx′(t,τ)对时延τ的一维傅氏变换。

将（5.5.39）式带入(5.5.38)式,得到

(5.5.42)(4)广义谱相关函数的对θ的一维傅里叶变换。

式中

(5.5.43)(5.5.44)Rx′(θ,f)称为广义谱相关函数,而

(5.5.45)Ψ(θ,f)称为谱相关域的核函数,由上式得到

(5.5.46)谱相关函数Rx(θ,ξ)为

(5.5.47)因此

(5.5.48)

上面四种等价的表示形式可归纳为对维格纳分布、模糊函数，瞬时相关函数和瞬间频相关函数的四种表示形式。它们是

(5.5.49)(5.5.50)(5.5.51)(5.5.52)四种核函数之间的关系为

(5.5.53)(5.5.54)(5.5.55)四种核函数之间的关系如图5.5.3所示。

图

5.5.2Cx(t,f)在(t,τ)、

(θ,f)、

(θ,τ)及(t,f)域之间的关系

图

5.5.3四种核函数之间的关系

WD是能量化的时频表示，而模糊函数是相关化的时频表示，相应的前者是时频域（t，f），后者是模糊域（θ,τ）。相应的可以将Cohen类时频分布分成两类时频表示，即能量化时频表示和相关化时频表示。能量化时频表示将瞬时功率和谱能量密度两种概念综合在一起，能量化的解释用边缘特性(5.5.21)式和(5.5.22)式表示，它的原始公式是(5.5.31)式，(5.5.49)～(5.5.52)式是它的四种等价形式。能量化时频表示简称CE类(下标E表示能量化)。相关化时频表示将瞬时相关和相关综合在一起，并用相关化边缘特性描述(即(5.5.26)、(5.5.27)式)，相关化时频表示简称Cc类(下标c表示相关化)。CE类时频表示Cx(t,f)中任何一种等价形式的对偶相关化时频表示Pdual,x(θ,τ)均有下面公式：

(5.5.56)CE类和Cc类的傅氏变换关系用图

5.5.4表示。

图

5.5.4CE类和Cc类的傅氏变换关系

5.5.3广义双线性时频分布举例

1.指数分布（ED）这种广义双线性时频分布的核函数是指数形式的，因此称为指数分布。又因为是由Choi

Williams提出的，也称Choi

Williams分布（CWD）。这种分布对含有多个频率分量的信号，可以有效地抑制交叉项，并保持时频分布的一些期望特性，是一种常用的广义双线性时频分布，下面作简单介绍。

指数分布在（θ,τ）域的核函数为

(5.5.57)在(t,τ)域的核函数为

(5.5.58)式中σ>0，是尺度因子。将上式带入（5.5.38）、(5.5.39)式中,得到它的时频分布Ex(t,f)为

(5.5.59)下面说明按照核函数的约束条件,这种时频分布具有的性质。(1)按照(5.5.57)式,核函数满足下式:因此该时频分布是实值的。

（2）

按照(5.5.57)式,核函数满足下列式:因此指数分布具有边缘特性和能量特性，

即

（3）信号的时移或频移在指数分布中产生相应的时移或频移，这是因为所有的CE类广义时频分布均有这种性质。

（4）指数分布对时间t的一阶矩和对频率的一阶矩分别等于群延迟T（Ω）和瞬时频率Ω（t）特性。这是因为核函数满足

（5）指数分布可以有效地抑制多频率成分信号的交叉项。假设信号是一个多频率信号，用下式表示:(5.5.60)可以推出它的广义模糊函数用下式表示:(5.5.61)上式中第一项是广义自模糊函数项，第二项是广义互模糊函数项，

它们分别用下式表示：l≠m(5.5.62)(5.5.63)广义自模糊函数主要集中在原点附近,而广义互模糊函数或者说交叉项则离开原点一段距离。这样,如果要突出广义自模糊函数项，并且抑制广义互模糊函数项，则要求核函数在（θ,τ）平面的原点附近有较大的加权值，对交叉项的加权值尽可能小，从而达到抑制交叉项的目的。指数分布的核函数则具备这种性质，因此它可以有效地抑制交叉项。假设两频率信号如下式所示：

(5.5.64)该信号的指数分布为

(5.5.64)式中加权值

(5.5.65)上面两式表明交叉项的幅度由WEIGHT控制,当Ω=(Ω1+Ω2)/2时,交叉项的幅度与成反比,另外,它以距离［Ω-(Ω1+Ω2)/2］2式指数地衰减。

2.广义指数分布（GED）与巴特沃斯分布（BUD）

1)广义指数分布广义指数分布的核函数为

（5.5.66）

式中,N、M为正的幂次;θ1和τ1分别是正的频率和时间尺度常数。下面说明该核函数具有低通特性。

先分析下面低通滤波函数θM(x)的特性:（5.5.67）

该函数的滤波特性如图5.5.5所示。此图表明,当x1一定时,随M的增大,通带变得愈平坦,过渡带愈窄。当M→∞时,变成在x=x1

处截止的理想低通滤波器。

图

5.5.5θM(x)的滤波特性

令 ,则

(5.5.68)这种情况下,广义指数核函数变成一维低通滤波函数θM(x),当N=M=1时, ,x=τθ,上式变为(5.5.69)即θ1(x)是指数分布的核函数。此时,核函数具有最差的低通特性。当M加大时，滤波性能变好。因此说,广义指数分布是可以改善时频分布性能的一种分布，通过改善,使其自项失真小，

交叉项衰减大。

2）巴特沃斯分布巴特沃斯分布的核函数定义为

N,M,θ1,τ1>0(5.5.70)上式表明,该核函数具有可变的、平坦的通带及窄的过渡带的低通特性。当 ,

上式变成

(5.5.71)上式即是滤波器中的最大平坦平方幅度的巴特沃斯滤波器的幅度平方函数,当M加大时，它趋近于理想低通滤波器。因此，这种分布也可以改善指数分布的时频分布的性能。除了上面介绍的广义双线性时频分布以外，常用的还有平滑伪维格纳分布（SPWD）、锥形（Cone）核分布（CKD）、降低干扰项分布（RID）、贝塞尔（Bessel）分布（BD）等等。

5.6时频分析在编队目标架次检测中的应用

5.6.1编队飞机目标多普勒特性分析本节我们分析编队飞机目标的多普勒频率特性，从而探讨在多普勒域实现多目标分辨的可能性。众所周知，雷达目标回波的多普勒频率fd取决于雷达工作波长λ、目标运动速度v以及目标飞行方向与雷达视线的夹角φ，

即

(5.6.1)

图5.6.1显示了两架飞机编队情况。假设编队飞行时，两架飞机的速度相同，则它们回波的多普勒频率的差别就取决于它们各自与雷达视线夹角的不同，

即

(5.6.2)φ1与φ2的差别是由两飞机的间距d和目标距雷达的距离R所决定的。一般有R>>d,

因此，φ1与φ2的差别是很小的，进一步化简，

有

(5.6.3)图5.6.1两架飞机编队情况式中，利用了当R>>d

时，sin(φ1-φ2)≈φ1-φ2≈(d/R)sinφ,其中φ=(φ1+φ2)/2是平均的雷达视线与飞行方向的夹角。由(5.6.3)式可以看出，