CN110471018B - 一种频谱校正方法 - Google Patents

Abstract

本发明利用变窗长方法，求解得到系列幅值、频率和相位。根据主瓣中心与最近谱线的距离随窗长变化的原理，选取不同窗长的窗函数进行截断，获得多组FFT变换结果，则最近谱线相对于主瓣中心的位置将会存在规律性的接近和远离。理论上，只要遍历各种窗函数长度，最近谱线与主瓣中心的距离就可以无限接近。此时，最大幅值对应的谱线位置即与主瓣中心最近，即最接近幅值真值，并且主瓣中心与最近谱线的距离为零均值的分布。因此，可以将系列频率值的均值作为校正的频率值，将系列幅值的最大值作为校正的幅值，将系列相位的均值作为校正的相位值。本发明频率值、幅值和相位值三者的校正过程相互独立，任一一者的校正误差，不会传递并影响其他两者。

Description

一种频谱校正方法

技术领域

本发明涉及数字信号处理领域，具体为一种频谱校正方法。

背景技术

离散频谱分析实现了信号处理从时域到频域的转变，推动了计算机应用技术的发展，在机械、电子、仪器仪表等领域得到了广泛应用。频谱分析准确性对工程应用具有十分重要的意义，然而计算机难以处理实际连续信号，需对信号进行截断和离散，当采样长度不是信号周期整数倍时，将会造成频谱泄漏。

频谱泄漏可分为长程谱泄漏和短程谱泄漏。从连续傅里叶变换到DFT，需要经过时域离散、数据截断以及频域离散过程。信号的时域离散化导致频域周期化，根据Nyquist采样定理，采样频率应大于信号最高频率的两倍，否则会产生假频造成频率混叠现象。计算机处理信号的长度总是有限的，必然要对信号进行截断，若截断信号长度N是信号周期的非整数倍，即非同步采样，在经频域离散造成时域周期化后，截断处会因吉布斯现象产生振荡，这种不连续状态将导致长程谱泄漏，泄漏程度与窗谱的旁瓣特性密切相关。在非同步采样中，频域离散还会使得信号真实频率f₀位于两条离散谱线k和k+1之间，造成短程谱泄漏，这种现象称为栅栏效应，相邻谱线间的宽度Δf＝f_s/N为频率分辨率，直接影响频谱分析精度。

这种频谱泄漏现象会影响频谱分析的准确性，给各类工程应用造成障碍。例如旋转机械振动响应信号包含了转频及其倍频成分，在机械故障诊断中需要研究各倍频轴心轨迹的特征，而短程谱泄漏造成的以DFT频谱分析为基础的频谱相位和幅值的不准确，影响了各频率成分的提取和提纯轴心轨迹的合成。

为减小短程谱泄漏造成的频谱分析误差，研究者提出了多种短程谱泄漏抑制方法，如频谱细化法、插值法、能量重心法、相位差法、三角形法等。上述方法能有效抑制频谱泄漏，进行相对精确的频谱分析和参数估计。但是频谱细化方法存在计算复杂度高或细化频带窄的缺点，而上述其他方法都基于主瓣中心与最近谱线的距离Δx进行频谱校正，Δx的计算误差将影响相位幅值和频率的校正，存在误差传递的问题。

发明内容

针对现有技术中存在的问题，本发明提供一种频谱校正方法，能够实现频谱幅值、频率和相位独立校正，不需要计算主瓣中心与最近谱线的距离Δx，避免了校正误差传递。

本发明是通过以下技术方案来实现：

一种频谱校正方法，包括如下步骤，

步骤1，对采样频率为f_s的离散正弦信号S(n)＝[s₀ s₁ s₂ … s_k]，以长度为N的窗函数W₀截取信号，获取第一短时信号y₀(n)＝[s₀ s₁ … s_N-1]；其中，k为信号数据点数，1<k<∞，1<N<k；

步骤2，通过FFT算法求解短时信号y₀(n)，得到其频域离散复数序列

通过

的最大绝对值序号j₀和频率分辨率Δf₀计算得到信号的第一近似频率值f_a0，并将序号j₀处

的绝对值记为第一近似幅值A_a0，将序号j₀处

的虚部与实部反正切值

记为第一相位值

步骤3，依次改变窗长为N+i，i＝1,2,…m，1<m<∞，重复步骤1和2获取m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

步骤4，将m个频率值的均值作为校正的频率值，将m个幅值的最大值最为校正的幅值，将m个相位的均值作为校正的相位值，完成对采样频率为f_s的离散正弦信号的频谱校正。

优选的，步骤2中，将最高谱线对应的频率值记为频率真值的近似值，即f_a0＝(j₀-1)Δf₀。

优选的，步骤2中，将最高谱线的幅值记为频率幅值的近似值，即

优选的，步骤2中，将最高谱线处复数相角记为相位值，即

优选的，步骤4中，校正的频率值为

优选的，步骤4中，校正的幅值为

优选的，步骤4中，校正相位值

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

本发明利用变窗长方法，求解得到系列幅值、频率和相位。根据主瓣中心与最近谱线的距离随窗长变化的原理，选取不同窗长的窗函数进行截断，获得多组FFT变换结果，则最近谱线相对于主瓣中心的位置将会存在规律性的接近和远离。理论上，只要遍历各种窗函数长度，最近谱线与主瓣中心的距离就可以无限接近。此时，最大幅值对应的谱线位置即与主瓣中心最近，即最接近幅值真值，并且主瓣中心与最近谱线的距离为零均值的分布。因此，可以将系列频率值的均值作为校正的频率值，将系列幅值的最大值作为校正的幅值，将系列相位的均值作为校正的相位值。本发明频率值f、幅值A和相位值

三者的校正过程相互独立，任一一者的校正误差，不会传递并影响其他两者。校正过程与窗函数的具体表达式无关，可以适应多种加窗信号。

附图说明

图1为本发明实例中所述的频谱校正方法流程图。

图2为本发明实例中所述的加窗获取短时信号示意图。

图3为本发明实例中所述的局部频带内最高谱线位置示意图。

图4为本发明实例中所述的实例信号时域波形。

图5为本发明实例中所述实例信号的短时信号频谱图。

图6为本发明实例中所述的第一频率成分幅值序列的分布图。

图7为本发明实例中所述的第一频率成分频率序列的分布图。

图8为本发明实例中所述的第一频率成分相位序列的分布图。

图9为本发明实例中所述的第二频率成分幅值序列的分布图。

图10为本发明实例中所述的第二频率成分频率序列的分布图。

图11为本发明实例中所述的第二频率成分相位序列的分布图。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明一种频谱校正方法，能够独立校正频谱幅值、频率和相位独立校正，包括如下步骤，

步骤1，对采样频率为f_s的离散正弦信号S(n)＝[s₀ s₁ s₂ … s_k]，(k为信号数据点数，1<k<∞)以长度为N(1<N<k)的窗函数W₀截取信号，获取第一短时信号y₀(n)＝[s₀ s₁ …s_N-1]；

步骤2，通过FFT算法求解短时信号y₀(n)，得到其频域离散复数序列

进而，通过

的最大绝对值序号j₀和频率分辨率Δf₀计算得到信号的第一近似频率值f_a0，并将j₀序号处

的绝对值记为第一近似幅值A_a0，将j₀序号处

的虚部与实部反正切值

记为第一相位值

步骤3，依次改变窗长为N+1,N+2,…,N+m，1<m<∞，可根据精度需求选取。采用步骤1和2获取m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

步骤4，将m个频率值的均值作为校正的频率值，将m个幅值的最大值最为校正的幅值，将m个相位的均值作为校正的相位值，完成对采样频率为f_s的离散正弦信号的频谱校正。

具体的，包括以下步骤：

1)以长度为N的窗函数W₀截取采样频率为f_s的离散正弦信号S(n)＝[s₀ s₁ s₂ …s_k]，获取短时信号y₀(n)＝[s₀ s₁ … s_N-1]；

2)对短时信号y₀(n)做FFT变换，得到其频域离散复数序列

3)寻找离散复数序列

绝对值的最大值所对应的序号，即最高谱线对应的序号，并将其记为j₀；

4)将最高谱线对应的频率值记为频率真值的近似值，即f_a0＝(j₀-1)Δf₀；将最高谱线的幅值记为频率幅值的近似值，即

将最高谱线处复数相角记为相位值，即；

5)改变窗长为N+i，i＝1,2,…m，重复步骤1)～4)，获得m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

6)将m个频率值的均值作为校正的频率值f_c，将m个幅值的最大值最为校正的幅值A，将m个相位的均值作为校正的相位值

其中，频率值f、幅值A和相位值

三者的校正过程相互独立，任一一者的校正误差，不会传递并影响其他两者。正过程与窗函数的具体表达式无关，可以适应多种加窗信号。

本发明所述的方法实施时的基本流程如图1所示，改变窗长为N+i，i＝1,2,…m，求解得到系列幅值、频率和相位，利用其分布特征，将系列频率值的均值作为校正的频率值，将系列幅值的最大值最为校正的幅值，将系列相位的均值作为校正的相位值。

首先，如图2所示，以长度为N的窗函数W₀截取采样频率为f_s的离散正弦信号S(n)＝[s₀ s₁ s₂ … s_k]，获取短时信号y₀(n)＝[s₀ s₁ … s_N-1]；

然后，对短时信号y₀(n)做FFT变换，得到其频域离散复数序列

再寻找离散复数序列

绝对值的最大值所对应的序号，即最高谱线对应的序号，并将其记为j₀，如图3所示；

将最高谱线对应的频率值记为频率真值的近似值，即f_a0＝(j₀-1)Δf₀；将最高谱线的幅值记为频率幅值的近似值，即

将最高谱线处复数相角记为相位值，即

然后，如图1所示，改变窗长为N+i，i＝1,2,…m，获得m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

即频率为f_ai＝(j₀-1)Δf₀，幅值为

相位为

最后，将m个频率值的均值作为校正的频率值f_c，将m个幅值的最大值作为校正的幅值A_c，将m个相位的均值作为校正的相位值

即校正频率为

校正幅值为A_c＝max(A_ai)，校正相位为

本优选实例，以信号

为例，其时域波形如图7所示。按照图1所示的基本流程，采样频率为f_s＝1000Hz，取窗长N＝1000，得到第一个短时信号的离散频谱如图8所示，按照上述实施方法，计算得到第一组两个频率成分的幅值、频率和相位值为A_a1＝4.92,12.85、f_a1＝53,79和

改变窗长为N+i，i＝1,2,…m，m＝1000分别获得两个频率成分的m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

其分布图分别如图6～图11所示。最后，将两个频率成分别分对应的m个频率值的均值作为校正的频率值f_c，将m个幅值的最大值最为校正的幅值A_c，将m个相位的均值作为校正的相位值

校正值、真值及误差如表1所示。为对比说明，给出了本例中直接FFT计算的最大误差值作为参照。从表1中可以看出，直接FFT变换中，频率成分1的幅值、频率和相位误差依次为0.7519、0.42、1.56749，频率成分2的幅值、频率和相位误差依次为2.2601、0.4047、1.56711。其中幅值和相位误差均较大，难以满足工程应用需求，频率误差较小。而校正后，频率成分1的幅值、频率和相位误差依次为-0.0049、-0.0367、-0.00205，频率成分2的幅值、频率和相位误差依次为-0.0141、-0.0538、-0.00423。校正后，幅值和相位误差显著降低，频率误差也进一步降低。

表1校正结果对比

Claims (1)

1.一种频谱校正方法，其特征在于，包括如下步骤，

步骤1，对采样频率为f_s的离散正弦信号S(n)＝[s₀ s₁ s₂ … s_k]，以长度为N的窗函数W₀截取信号，获取第一短时信号y₀(n)＝[s₀ s₁ … s_N-1]；其中，k为信号数据点数，1<k<∞，1<N<k；

步骤2，通过FFT算法求解短时信号y₀(n)，得到其频域离散复数序列

通过

的最大绝对值序号j₀和频率分辨率Δf₀计算得到信号的第一近似频率值f_a0，并将序号j₀处

的绝对值记为第一近似幅值A_a0，将序号j₀处

的虚部与实部反正切值

记为第一相位值

步骤3，依次改变窗长为N+i，i＝1,2,…m，1<m<∞，重复步骤1和2获取m个频率值f_a1,f_a2，…，f_am、幅值A_a1,A_a2，…，A_am和相位值

步骤4，将m个频率值的均值作为校正的频率值，将m个幅值的最大值最为校正的幅值，将m个相位的均值作为校正的相位值，完成对采样频率为f_s的离散正弦信号的频谱校正；

校正的频率值为

校正的幅值为A_c＝max(A_ai)；

校正相位值

Images