CN110060700A - 基于参数谱估计的短序列音频分析方法 - Google Patents

Abstract

本发明公开了一种基于参数谱估计的短序列音频分析方法，包括：建立音频信号x(n)的数学表达公式；利用信号构造Hankel矩阵Y；计算Y的协方差矩阵R，并将R做奇异值分解得到酉矩阵U、V和奇异值矩阵∑；利用矩阵∑中奇异值计算信号的奇异熵增量，从而判别信号模态阶数g；根据模态阶数将矩阵V划分成信号空间V₁和噪声空间V₂，并从信号空间中取出两个子矩阵V₃和V₄；求解V₃和V₄间的映射关系ψ；根据矩阵ψ计算信号x(n)的频率分量和对应的幅度、相位；估计信号噪声功率；计算信号信纳比、失真度。本发明结合基于奇异熵的信号模态阶数识别算法分离信号空间，计算各模态的频率、幅度、相位信息和信号噪声功率，而计算信纳比、失真度，测量精度高。

Description

基于参数谱估计的短序列音频分析方法

技术领域

本发明属于通信测试领域，尤其涉及到一种基于参数谱估计的短序列音频分析方法，可用于便携式无线电综合测试仪对音频信号的测量分析。

背景技术

在通信测试中，通常采样导频频率、导频频偏、发音频谐波失真等指标衡量发射机性能，采用模拟灵敏度、收音频谐波失真等指标衡量接收机性能。上述指标的测量过程分别依赖于对音频信号频率、功率、信纳比、失真度等参数的测量，即音频分析。在传统测试系统中，多采用无线电综合测试仪的音频分析功能实现这部分测试需求，然而随着通信测试领域的不断细分，出现的一些现场级测试以及快速维修保障需求在设备的便携性和经济性方面对测试仪器提出了更高要求，此时笨重而昂贵的无线电综合测试仪不在成为首选。因此越来越多的方案更倾向于采用信号采集卡配合便携式计算机搭建测试系统，而音频分析工作则全权交由软件负责。

音频分析技术是一种通过波形分析方法测量音频信号各正弦分量频率、功率、相位以及电平值、信纳比、失真度、信噪比等特征的技术，目前大多采用基于FFT的方法实现。该类方法简单且计算代价较小，够能实现对整周期单音信号的高精度分析。然而由于离散傅里叶变换存在显著的频谱泄露和栅栏效应，该类方法在测量存在多次间谐波的复杂信号时精度大幅下降。另一方面，基于FFT的参数估计方法受到不确定原理的傅里叶分辨率的限制，而实际应用场景中常出现频率靠近基波的谐波分量，受限于常用音频采集卡的采样长度，基于FFT的方法往往无法分辨这类信号。

发明内容

本发明的目的是克服已有音频分析方法中测量精度和频率分辨率受限于采样序列长度的缺点，突破普通音频采集卡连续采样时长的限制，提出一种基于参数谱估计的短序列音频分析方法。

为实现上述目的，本发明根据如下所述技术方案实现：

基于参数谱估计的短序列音频分析方法，其特征在于，包含以下步骤：

步骤S1：信号采集卡输入或软解调器输出音频信号有如下表示形式：

其中K为信号总正弦分量数，第i次谐波的幅度、频率、相位分别由a_i、ω_i，表示，e(n)表示零均值的高斯白噪声；

步骤S2：引入1×m维向量Y′，并由X(n)构造Hankel矩阵Y：

Y′(n)＝[x(n),x(n+1),…,x(n+m-1)]

Y＝[Y′(0),Y′(1),…,Y′(L)]^T

其中M＝N-L；

步骤S3：计算Y的协方差矩阵R，并将方阵R^(N-L)×^(L+1)做奇异值分解：

R＝YY^T

R＝U∑V^T

其中对角阵∑为方阵R的奇异值矩阵，其对角元素d_ii为R的第i个奇异值，其中1≤i≤h,h＝min(N-L,L+1)；

步骤S4：计算信号的模态阶数，由奇异值矩阵∑计算奇异熵ΔE_i：

以ΔE_i拐点对应的阶数g作为信号的模态阶数，该拐点的识别方法为计算ΔE_i的微分ΔE_i′，并将满足ΔE_i′≥Threshould的最大i值作为信号的模态阶数；

步骤S5：根据确定的模态阶数g将矩阵V划分成两个子空间：

V＝[V₁ V₂]

其中V₁ ^M×2g为信号子空间，V₂ ^M×(M-2g)为噪声子空间，利用下面关系式在V₁中找出两个旋转子矩阵V₃和V₄：

V₃＝[I_M-10₁]·V₁

V₄＝[0₁I_M-1]·V₁

其中I_M-1为(M-1)×(M-1)的单位矩阵，0₁为(M-1)×1全零列向量；

步骤S6：根据移不变性质，有：

V₄＝V₃·ψ

求解ψ矩阵：

ψ＝(V₃ ^HV₃)^-1V₃ ^HV₄；

步骤S7：计算输入信号的频率分量以及对应的幅度和相位；

步骤S8：根据各模态的频率、幅度、相位值拟合原始输入信号x(n)，得到：

计算x′(n)和x(n)的均方误差e′²作为对噪声e(n)功率估计；

步骤S9：计算输入信号信纳比为：

计算输入信号失真度为：

其中a′_max＝max{a′_i|i＝1,2,…,P}。

上述技术方案中，通过输入信号构造矩阵Y的方式。

上述技术方案中，计算奇异熵并将奇异熵拐点的阶数作为对输入信号的模态阶数，模态阶数判别时的阈值决定了信号子空间和噪声子空间的分割。

上述技术方案中，利用已经求得的各模态频率、幅度、相位重构输入信号x(n)，并将重构信号和原始信号间的均方误差值作为对噪声功率的估计。

上述技术方案中，步骤S7中，信号的频率分量由下式计算：

其中λ_k为矩阵ψ的第p个特征值，f_p正负成对出现，剔除负频率后即所求信号频率分量，

各频率分量对应的幅度和相位计算方法如下，首先构造矩阵Λ：

Λ＝[λ_ψ ⁰,λ_ψ ¹,…,λ_ψ ^N-1]^T

其中λ_ψ ^m＝[λ₁ ^m,λ₂ ^m,…,λ_P ^m]，λ₁,λ₂,…,λ_P为矩阵ψ的P个特征值，复数域模态系数A_k满足:

ΛA＝X

其中1≤k≤P，A＝[A₁,A₂,…,A_P]^T,X＝[x(0),x(1),…,x(N-1)]^T，各频率分量的幅度和相位由下式计算：

a′_p＝2|A_p|,p＝1,2,…,P

与现有技术相比，本发明具有如下的有益效果：

本发明提出的基于参数谱估计的短序列音频分析方法通过算法推导和工程实现验证过的，可以有效分析较短的音频数据片段，计算其基波和各次谐波频率、功率、相位以及音频信纳比、失真度等信息，且具有较高的分析精度。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为本发明中音频信号的时域波形图；

图2为本发明的音频分析方法流程示意图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。

在模拟通信系统中，调制信号经过信道传输和调制解调后，往往会引入多种噪声与失真，而这些噪声与失真最终将表现在解调输出的音频信号上。因此通信测试时，常见的做法是由测试设备发生一单音调制信号输入待测接收机，通过观测其解调输出的音频信号基波频率、基波功率、信纳比、失真度等参数来衡量待测接收机性能指标；或者反向地，通过待测发射机发生一单音调制信号输入测试设备，通过观测测试设备解调输出的音频信号基波频率、基波功率、信纳比、失真度等参数来衡量待测发射机性能指标。当然，上述测试用例仅仅是本发明的一部分实施例，而不是全部实施例。

图1为本发明的音频信号的时域波形图，图2为本发明的音频分析流程图，如图2所示，本发明的基于参数谱估计的短序列音频分析方法，包含以下步骤：

步骤S1：信号采集卡输入或软解调器输出音频信号有如下表示形式：

其中K为信号总正弦分量数，第i次谐波的幅度、频率、相位分别由a_i、ω_i，表示，e(n)表示零均值的高斯白噪声；

步骤S2：引入1×m维向量Y′，并由X(n)构造Hankel矩阵Y：

Y′(n)＝[x(n),x(n+1),…,x(n+m-1)]

Y＝[Y′(0),Y′(1),…,Y′(L)]^T

其中M＝N-L；

步骤S3：计算Y的协方差矩阵R，并将方阵R^(N-L)×^(L+1)做奇异值分解：

R＝YY^T

R＝U∑V^T

其中对角阵∑为方阵R的奇异值矩阵，其对角元素d_ii为R的第i个奇异值，其中1≤i≤h,h＝min(N-L,L+1)；

步骤S4：计算信号的模态阶数，由奇异值矩阵∑计算奇异熵ΔE_i：

以ΔE_i拐点对应的阶数g作为信号的模态阶数，该拐点的识别方法为计算ΔE_i的微分ΔE_i′，并将满足ΔE_i′≥Threshould的最大i值作为信号的模态阶数；

步骤S5：根据确定的模态阶数g将矩阵V划分成两个子空间：

V＝[V₁ V₂]

其中V₁ ^M×2g为信号子空间，V₂ ^M×(M-2g)为噪声子空间，利用下面关系式在V₁中找出两个旋转子矩阵V₃和V₄：

V₃＝[I_M-10₁]·V₁

V₄＝[0₁I_M-1]·V₁

其中I_M-1为(M-1)×(M-1)的单位矩阵，0₁为(M-1)×1全零列向量；

步骤S6：根据移不变性质，有：

V₄＝V₃·ψ

求解ψ矩阵：

ψ＝(V₃ ^HV₃)^-1V₃ ^HV₄；

步骤S7：计算输入信号的频率分量以及对应的幅度和相位；

步骤S8：根据各模态的频率、幅度、相位值拟合原始输入信号x(n)，得到：

计算x′(n)和x(n)的均方误差e′²作为对噪声e(n)功率估计；信号的频率分量由下式计算：

其中λ_k为矩阵ψ的第p个特征值，f_p正负成对出现，剔除负频率后即所求信号频率分量，

各频率分量对应的幅度和相位计算方法如下，首先构造矩阵Λ：

Λ＝[λ_ψ ⁰,λ_ψ ¹,…,λ_ψ ^N-1]^T

其中λ_ψ ^m＝[λ₁ ^m,λ₂ ^m,…,λ_P ^m]，λ₁,λ₂,…,λ_P为矩阵ψ的P个特征值，复数域模态系数A_k满足:

ΛA＝X

其中1≤k≤P，A＝[A₁,A₂,…,A_P]^T,X＝[x(0),x(1),…,x(N-1)]^T，各频率分量的幅度和相位由下式计算：

a′_p＝2|A_p|,p＝1,2,…,P

步骤S9：计算输入信号信纳比为：

计算输入信号失真度为：

其中a′_max＝max{a′_i|i＝1,2,…,P}。

上述技术方案中，通过输入信号构造矩阵Y的方式，L的取值范围在N/4到N/3之间。

上述技术方案中，计算奇异熵并将奇异熵拐点的阶数作为对输入信号的模态阶数，模态阶数判别时的阈值决定了信号子空间和噪声子空间的分割。

上述技术方案中，利用已经求得的各模态频率、幅度、相位重构输入信号x(n)，并将重构信号和原始信号间的均方误差值作为对噪声功率的估计。

本发明基于参数谱估计的算法思想，结合基于奇异熵的信号模态阶数识别算法分离信号空间，计算各模态的频率、幅度、相位信息和信号噪声功率，从而计算信纳比、失真度，测量精度高，无分辨率限制，可适用于短序列音频信号分析。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。

Claims (5)

1.基于参数谱估计的短序列音频分析方法，其特征在于，包含以下步骤：

步骤S1：信号采集卡输入或软解调器输出音频信号有如下表示形式：

其中K为信号总正弦分量数，第i次谐波的幅度、频率、相位分别由a_i、ω_i，表示，e(n)表示零均值的高斯白噪声；

步骤S2：引入1×m维向量Y′，并由X(n)构造Hankel矩阵Y：

Y′(n)＝[x(n),x(n+1),…,x(n+m-1)]

Y＝[Y′(0),Y′(1),…,Y′(L)]^T

其中M＝N-L；

步骤S3：计算Y的协方差矩阵R，并将方阵R^(N-L)×(L+1)做奇异值分解：

R＝YY^T

R＝U∑V^T

其中对角阵∑为方阵R的奇异值矩阵，其对角元素d_ii为R的第i个奇异值，其中1≤i≤h,h＝min(N-L,L+1)；

步骤S4：计算信号的模态阶数，由奇异值矩阵∑计算奇异熵ΔE_i：

以ΔE_i拐点对应的阶数g作为信号的模态阶数，该拐点的识别方法为计算ΔE_i的微分ΔE_i′，并将满足ΔE_i′≥Threshould的最大i值作为信号的模态阶数；

步骤S5：根据确定的模态阶数g将矩阵V划分成两个子空间：

V＝[V₁ V₂]

其中V₁ ^M×2g为信号子空间，V₂ ^M×(M-2g)为噪声子空间，利用下面关系式在V₁中找出两个旋转子矩阵V₃和V₄：

V₃＝[I_M-1O₁]·V₁

V₄＝[O₁I_M-1]·V₁

其中I_M-1为(M-1)×(M-1)的单位矩阵，O₁为(M-1)×1全零列向量；

步骤S6：根据移不变性质，有：

V₄＝V₃·ψ

求解ψ矩阵：

ψ＝(V₃ ^HV₃)^-1V₃ ^HV₄；

步骤S7：计算输入信号的频率分量以及对应的幅度和相位；

步骤S8：根据各模态的频率、幅度、相位值拟合原始输入信号x(n)，得到：

计算x′(n)和x(n)的均方误差e′²作为对噪声e(n)功率估计；

步骤S9：计算输入信号信纳比为：

计算输入信号失真度为：

其中a′_max＝max{a′_i|i＝1,2,…,P}。

2.根据权利要求1所述的基于参数谱估计的短序列音频分析方法，其特征在于，步骤S2中，通过输入信号构造矩阵Y的方式。

3.根据权利要求1所述的基于参数谱估计的短序列音频分析方法，其特征在于，步骤S4中，计算奇异熵并将奇异熵拐点的阶数作为对输入信号的模态阶数，模态阶数判别时的阈值决定了信号子空间和噪声子空间的分割。

4.根据权利要求1所述的基于参数谱估计的短序列音频分析方法，其特征在于，步骤S8中，利用已经求得的各模态频率、幅度、相位重构输入信号x(n)，并将重构信号和原始信号间的均方误差值作为对噪声功率的估计。

5.根据权利要求1所述的基于参数谱估计的短序列音频分析方法，其特征在于，步骤S7中，信号的频率分量由下式计算：

其中λ_k为矩阵ψ的第p个特征值，f_p正负成对出现，剔除负频率后即所求信号频率分量，

各频率分量对应的幅度和相位计算方法如下，首先构造矩阵Λ：

Λ＝[λ_ψ ⁰,λ_ψ ¹,…,λ_ψ ^N-1]^T

其中λ_ψ ^m＝[λ₁ ^m,λ₂ ^m,…,λ_P ^m]，λ₁,λ₂,…,λ_P为矩阵ψ的P个特征值，复数域模态系数A_k满足:

ΛA＝X

其中1≤k≤P，A＝[A₁,A₂,…,A_P]^T,X＝[x(0),x(1),…,x(N-1)]^T，各频率分量的幅度和相位由下式计算：

a′_p＝2|A_p|,p＝1,2,…,P