CN102269803B - 基于时间延迟的离散频谱低频成分的校正方法 - Google Patents

Abstract

本发明涉及一种基于时间延迟的离散频谱低频成分的校正方法，其特征在于包括以下步骤：(1).将离散时间信号按照一定时间延迟分成等长度的三段，对每段时间信号分别求和，并获得三个求和值；(2).利用离散时间信号频谱的对称性，并利用傅立叶变换公式获得离散频谱低频成分的有关频谱校正的表达形式；(3).将三个求和值结合离散频谱低频成分的有关频谱校正的表达形式获得离散频谱中低频成分的频率和相位；根据所获得的频率和相位结合(2)频谱校正的表达形式并获得低频成分的幅值；(4).对相位进行修正，获得相位的校正值。其计算量小，计算精度高，简单易行，适用于任何窗函数，提高了离散频谱低频成分校正的效率和准确性。

Description

基于时间延迟的离散频谱低频成分的校正方法

技术领域

本发明涉及对离散频谱低频成分校正方法，特别地，涉及一种基于时间延迟的离散频谱低频成分的校正方法。

背景技术

在频谱校正领域，信号频率的高低是相对而言的，当频率小于频谱的两倍频率分辨率时，视为低频信号，反之视为高频信号。例如，对于30Hz单频信号，若频率分辨率为20Hz，则该信号在频谱校正领域应该视为低频信号；若频率分辨率为1Hz，该信号在频谱校正领域应该视为高频信号。

在工程实践中，对于低频信号，采集到的样本内没有足够多的波动周期数，频谱会受到负频率的干涉，造成频谱失真。在机械设备故障诊断领域，由涡流传感器检测到的大型旋转机械的轴振动信号一般为周期信号和准周期信号，且具有较高的信噪比，通常采用快速傅立叶变换来进行频谱分析。

低速重载设备如炼铁厂高炉炉顶布料齿轮箱、炼钢厂的转炉耳轴、大型化工反应釜等都有必要的低频测量。在生产中有些设备间歇是运行的，运转时间可能比所需的采样时间短，采集到的周期数不够。

在医疗、生物和化工等工业领域，经常需要用到高精度的相位差流量计，为了提高处理速度，或者为了实时检测流量计的转速和相位差信息，则希望数据采集时间尽可能的短，然而实际中常常流量计本身的转速较小，导致采样时间内采集的流量计旋转周期数不足。

由于采样周期数不足，信号频率可能小于两倍频率分辨率，频谱会受到负频率成分的干涉，使谱线的幅值和相位产生很大误差。

从理论上分析，对于单谐波高频率信号，加矩形窗时幅值的最大误差可达36.4％，即使加其它窗函数，也不能完全消除此误差，如加汉宁窗时，幅值的最大误差仍高达15.3％，相位误差高达90°。而对于低频信号，由于频谱受到了负频率的干涉，误差更大，校正更困难。

传统的频谱校正方法基本上均未考虑频谱中负频率的影响，当振动信号的频率很低时，负频率对频谱的干涉现象非常明显，使得传统校正算法精度明显下降，甚至失效。

国外学者没有对频谱校正开展过研究，国内对频谱校正的研究很多，但是对低频成分校正问题进行研究的学者非常少。根据文献《振动工程学报》2008.21(1)：第38-42页，记录文章《低频成分的频谱校正》作者陈奎孚、王建立、张森文提出了一种只适用于矩形窗的校正方法。它用的是频谱峰值附近的三条谱线，由于在大部分情况下三条谱线中有一条或两条幅值较小，所以该方法易受噪声的影响。根据文献《计量学报》2008.29(2)：第168-171页，记录文章《基于FFT的一种计及负频率影响的相位差测量新方法》作者张海涛、涂亚庆提出了基于相位差的低频校正方法，但校正算法较复杂，需要求解方程组，计算量较大。

本发明提出了一种区别于以上两种方法的离散频谱低频成分校正方法，通过对离散时间信号按照一定时间延迟进行适当的分段，对信号中低频成分的频率、相位和幅值信息进行精确的校正。

发明内容

本发明的目的在于，提出了一种离散频谱低频成分校正方法，通过对离散时间信号按照一定时间延迟进行适当的分段，实现了对离散时间信号中低频成分的频率、相位和幅值信息进行精确的地校正。其计算量小，计算精度高，简单易行，适用于任何窗函数，提高了离散频谱低频成分校正的效率和准确性。具体技术方案如下：

一种基于时间延迟的离散频谱低频成分的校正方法，包括以下步骤：

(1).将离散时间信号按照一定时间延迟分成等长度的三段，对每段时间信号分别求和，并获得三个求和值；

(2).利用离散时间信号频谱的对称性，使用傅立叶变换公式进行频谱分析，以便事先得到离散频谱低频成分的有关频谱校正的表达形式；

(3).将(1)所述的三个求和值结合(2)所述的离散频谱低频成分的有关频谱校正的表达形式中涉及低频成分的频率和相位的形式，获得离散频谱中低频成分的频率和相位，根据所述的离散频谱中低频成分的频率和相位结合(2)所述的离散频谱低频成分的有关频谱校正的表达形式中涉及低频成分的幅度的表达形式，并获得低频成分的幅值；

(4).根据(3)所述的频率、幅值和相位结合(2)所述的有关频谱校正的表达形式，对相位进行修正，获得相位的校正值。

获得所述的步骤(1)所述的三个求和值，包括以下步骤：

(1.1)将待分析的离散时间信号设为x(t)，采样长度为2N，采样时间为2T；

(1.2)截取该信号[T/2，3T/2)区间上的部分作为x₁(t)，[0，T)区间上的部分作为x₂(t)，[T，2T)区间上的部分作为x₃(t)；

(1.3)则每段信号时长为T，长度为N，相邻信号的时间延迟为T/2；

(1.4)对每段时间信号分别求和，获得三个求和值：其中，X_i，i＝1，2，3表示三段时间信号的求和值。

获得步骤(2)所述的离散频谱低频成分的有关频谱校正的表达形式的方法如下：

(2.1).设三段信号的傅立叶变换为X_i(f)，i＝1，2，3，f表示频率；

(2.2).利用离散时间信号频谱的对称性，使用傅立叶变换公式，可获得离散频谱低频成分的频谱校正的表达形式为：

$\cos \left(\mathrm{\π}{\mathrm{Tf}}_0\right)=\frac{X_2\left(0\right)+X_3\left(0\right)}{2X_1\left(0\right)}$

$\tan \mathrm{\θ}=\frac{X_3\left(0\right)\cos \left(\mathrm{\πT}{f}_0\right)-X_1\left(0\right)}{X_3\left(0\right)\sin \left(\mathrm{\π}{\mathrm{Tf}}_0\right)}$

$A=\frac{X_3\left(0\right)}{W\left(f_0\right)\cos \left(\mathrm{\θ}\right)}$

式中，f₀表示离散时间信号中低频成分的频率校正值，X_i(0)，i＝1，2，3分别表示三段信号频谱0Hz处的谱线值，θ表示低频成分的相位校正值，A表示低频成分的幅值校正值，W(f₀)表示窗函数频谱在频率f₀处的幅值；

(2.3).利用离散傅立叶变换公式的性质，可获得X_i(0)与X_i，i＝1，2，3之间的关系为：X_i(0)＝X_i，i＝1，2，3.。

获得步骤(3)所述的离散频谱中低频成分的幅值、频率和相位的校正值的具体步骤如下：

根据X_i(0)和X_i，i＝1，2，3相等，可将X_i，i＝1，2，3结合频谱校正的表达形式，得到频率f₀、相位θ和幅值A，其中频率f₀和幅值A为所要求的校正值，相位θ经过修正公式后得到所要求的相位校正值。

所述的修正公式为θ′＝θ+180f₀/f_s，其中，f_s为采样频率，θ′即为所要求的相位校正值。

所述的本方法适用于对称窗函数。

综上所述，本方法简单易行，适用于任何对称窗函数，对离散频谱上小于两个频率分辨率的低频频成分的幅值、频率和相位，有较高的校正精度。有效弥补了已有方法不能校正低频成分的不足。本发明具有下列区别于其它方法的显著优势：

1)使得获得的频谱校正公式简单、参数少和鲁棒性强；

2)不需要进行快速傅立叶变换或其它变换，计算量小；

3)本发明提出的方法适用于任何对称窗函数，适用范围广；

4)本发明简单可靠，精度高，便于工程实践中使用。

附图说明

图1是本发明涉及实施例的信号分段简图。

图2是本发明涉及实施例的相位误差与频率关系图。

图3是本发明涉及实施例的校正误差与频率关系图。

图4是本发明涉及实施例的校正误差与相位关系图一。

图5是本发明涉及实施例的校正误差与相位关系图二。

具体实施方式

下面结合实施例进一步描述本发明。本发明的范围不受这些实施例限制。结合附图对本发明的内容作进一步详细说明如下：

1)将离散时间信号按照一定时间延迟分成等长度的三段，对每段时间信号分别求和，获得三个求和值；

2)利用离散时间信号频谱的对称性，从傅立叶变换公式入手，推导出离散频谱低频成分的频谱校正公式；

3)将三个求和值带入频谱校正公式，获得离散频谱中低频成分的幅值、频率和相位的校正值。

时窗函数一般为实偶函数，即w(t)＝w(-t)，其傅氏变换W(f)也为实偶函数。对信号x(t)加长度为T的窗函数w_T(t)得y(t)＝x(t)w_T(t)。其中，w_T(t)由窗函数w(t)在时间上平移T/2得到，即w_T(t)＝w(t-T/2)。设窗函数w(t)的傅立叶变换为F[w(t)]＝W(f)，由傅立叶变换的时移特性可得F[w_T(t)]＝W(f)e^-jπfT。设有一谐波信号x(t)＝Acos(2πf₀t+θ)，A，f₀和θ分别为幅值，频率和初相位。采样长度为2N，采样时间为2T。参照图1所示，截取该信号[T/2，3T/2)区间上的部分作为x₁(t)，[0，T)区间上的部分作为x₂(t)，[T，2T)区间上的部分作为x₃(t)。则每段信号时长为T，长度为N，相邻信号的时间延迟为T/2。

将x₁(t)乘以窗函数后，做傅立叶变换获得：

$X_1\left(f\right)=\frac{A}{2}W(f+f_0)e^{-j[\mathrm{\πT}(f+f_0)+\mathrm{\θ}]}+\frac{A}{2}W(f-f_0)e^{-j[\mathrm{\πT}(f-f_0)-\mathrm{\θ}]}---\left(1\right)$

将x₂(t)乘以同样的窗函数获得：

$X_2\left(f\right)=\frac{A}{2}W(f+f_0)e^{-j[\mathrm{\πTf}+2\mathrm{\π}{f}_{0}T+\mathrm{\θ}]}+\frac{A}{2}W(f-f_0)e^{-j[\mathrm{\πTf}-2\mathrm{\π}{f}_{0}T-\mathrm{\θ}]}---\left(2\right)$

将x₃(t)乘以同样的窗函数获得：

$X_3\left(f\right)=\frac{A}{2}W(f+f_0)e^{-j[\mathrm{\πTf}+\mathrm{\θ}]}+\frac{A}{2}W(f-f_0)e^{-j[\mathrm{\πTf}-\mathrm{\θ}]}---\left(3\right)$

这里通过时间延迟从x₃(t)相位中消去f₀将简化后面的计算。对于低频信号，在非整周期采样的情况下，频谱由于存在泄漏必然在0Hz处有值，所以将f＝0代入以式(1)-(3)获得：

$X_1\left(0\right)=\frac{A}{2}W\left(f_0\right)e^{-j[\mathrm{\π}{\mathrm{Tf}}_0+\mathrm{\θ}]}+\frac{A}{2}W\left({-f}_0\right)e^{-j[\mathrm{\πT}\left({-f}_0\right)-\mathrm{\θ}]}=\mathrm{AW}\left(f_0\right)\cos (\mathrm{\π}{\mathrm{Tf}}_0+\mathrm{\θ})---\left(4\right)$

$X_2\left(0\right)=\frac{A}{2}W\left(f_0\right)e^{-j[\mathrm{\πT}{f}_0+\mathrm{\θ}+\mathrm{\π}{f}_{0}T]}+\frac{A}{2}W\left({-f}_0\right)e^{-j[\mathrm{\πT}\left({-f}_0\right)-\mathrm{\θ}-\mathrm{\π}{f}_{0}T]}=\mathrm{AW}\left(f_0\right)\cos (2\mathrm{\π}{\mathrm{Tf}}_0+\mathrm{\θ})---\left(5\right)$

$X_3\left(0\right)=\frac{A}{2}W\left(f_0\right)e^{-j[\mathrm{\πT}{f}_0+\mathrm{\θ}-\mathrm{\π}{f}_{0}T]}+\frac{A}{2}W\left({-f}_0\right)e^{-j[\mathrm{\πT}\left({-f}_0\right)-\mathrm{\θ}+\mathrm{\π}{f}_{0}T]}=\mathrm{AW}\left(f_0\right)\cos \left(\mathrm{\θ}\right)---\left(6\right)$

这里应用了傅立叶变换的奇偶性质，当w(t)是实偶函数时，W(f)也为实偶函数。

由(4)-(6)式推导可以获得：

$\cos \left(\mathrm{\πT}{f}_0\right)=\frac{X_2\left(0\right)+X_3\left(0\right)}{2X_1\left(0\right)}---\left(7\right)$

$\tan \mathrm{\θ}=\frac{X_3\left(0\right)\cos \left(\mathrm{\πT}{f}_0\right)-X_1\left(0\right)}{X_3\left(0\right)\sin \left(\mathrm{\πT}{f}_0\right)}---\left(8\right)$

这里(7)，(8)式即为离散频谱中低频成分的频率和相位校正公式，求得f₀和θ后，带入(6)式即可获得幅值A

$A=\frac{X_3\left(0\right)}{W\left(f_0\right)\cos \left(\mathrm{\θ}\right)}---\left(9\right)$

其中X₁(0)，X₂(0)和X₃(0)分别为x₁(t)，x₂(t)和x₃(t)作离散傅立叶变换后的第一条谱线。一个长为N点的信号x(t)的离散傅立叶变换公式为：

$X\left(l\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_2\left(n\right)e^{-j2\mathrm{\π}\ln /N},l=\mathrm{0,1}...,N-1---\left(10\right)$

式中，l表示离散频谱的谱线号。代入l＝0，可获得X_i(0)，i＝1，2，3(i＝1，2，3)：

$X_i\left(0\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(n\right)\exp (-j2\·0\·n/N)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(n\right),i=\mathrm{1,2,3}---\left(11\right)$

用X_i，i＝1，2，3表示三段时间信号的求和值，即

$X_i=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(n\right),i=\mathrm{1,2,3}---\left(12\right)$

根据统计研究，发现由(8)式所获得的相位θ与其理论值之间存在一个有规律变化的偏差，该偏差仅与f₀与采样频率f_s有关，与幅值和相位无关，是关于归一化频率的过零点的线性函数，如图2所展示。图2中，纵轴为相位误差大小，横轴为归一化频率，斜率为k＝-180。为消除该偏差，应在用(7)-(9)式计算出频率、幅值和相位后对相位θ进行修正，才能获得理想的相位校正值：

θ′＝θ-kf₀/f_s                 (13)

以上所有推导没有具体利用哪一种窗函数，所以该方法适用于所有的对称窗函数。

以上方法可简单概述为，将一段时长为2T的序列，分成3段长为T，时差为T/2的序列后分别求和，结果带入式(7)-(9)求解频率、相位和幅值，最后修正相位。式(8)可进一步变成

$\tan \mathrm{\θ}=\cot \left(\mathrm{\πT}{f}_0\right)-\frac{X_1\left(0\right)}{X_3\left(0\right)\sin \left(\mathrm{\πT}{f}_0\right)}---\left(14\right)$

由于式(14)其中含有cot(πTf₀)，而余切函数为周期函数，周期为π，在0和π处为无穷，因而这里πTf₀的范围为0＜πTf₀＜π，即0＜f₀＜1/T＝f_s/N。所以f₀的取值范围大于0小于一个信号分段的频率分辨率。注意，这里的N和T为三段序列的长度和时间。所以f₀的取值范围实际上是小于两倍原始信号x(t)的频率分辨率(f_s/N＝2(f_s/2N))。

三段信号的时间关系图如图1所示。由于(8)式求得的相位为x₁(t)的相位，从图中可以看出，最后求得的修正相位为整个信号x(t)的1/4长度处的相位。

本实施例应用于验证基于时间延迟的离散频谱低频成分校正方法的精度。用计算机仿真信号x(t)＝cos(2πft+θ)，采样频率为2048Hz，采样点数1024，频率分辨率为2Hz，选用矩形窗。从信号x(t)的第513个点截取到第1536个点，作为x₁(t)；第1个点截取到第1024个点，作为x₂(t)；第1025个点截取到第2048个点，作为x₃(t)。利用(12)式对三段信号分别求和，获得X_i，i＝1，2，3。因为X_i＝X_i(0)，i＝1，2，3，所以可将X_i代入(7)-(9)式求解频率f₀、相位θ和幅值A。获得的频率f₀和幅值A为所求的校正值。相位θ还需要通过式(13)进行修正，获得θ′，θ′即为相位校正值。

图3展示的是x(t)中的f变化时，该方法对频率f₀、相位θ和幅值A的校正误差与f的关系。如图2所示，f从0.2Hz扫描到3.8Hz，步长为0.08Hz，初相位θ为60°。图4和图5展示的是x(t)中的θ变化时，该方法对频率f₀、相位θ和幅值A的校正误差与θ的关系。如图4和图5所示，x(t)的频率为3Hz，相位从0°扫描到179°，步长为1°。这里只考虑无噪声干扰的情形。

参照图3-图5，我们可以看出：新方法对低频信号的频率和相位校正精度非常高，接近双精度运算的下限，且误差基本不随频率和相位的变化而变化；幅值的校正误差相对稍高，误差不随相位的变化而变化，但随频率的增大而增大，但最大处误差依然小于0.001％。

Claims (4)

1.一种基于时间延迟的离散频谱低频成分的校正方法，其特征在于包括以下步骤：

(1).将离散时间信号按照一定时间延迟分成等长度的三段，对每段时间信号分别求和，并获得三个求和值；

(2).利用离散时间信号频谱的对称性，使用傅立叶变换公式进行频谱分析，以便事先得到离散频谱低频成分的有关频谱校正的表达形式；

(3).将(1)所述的三个求和值结合(2)所述的离散频谱低频成分的有关频谱校正的表达形式中涉及低频成分的频率和相位的表达式，获得离散频谱中低频成分的频率和相位，再根据所得到频率和相位结合(2)所述的离散频谱低频成分的有关频谱校正的表达形式中涉及低频成分的幅度的表达形式，并获得低频成分的幅值；

(4).根据(3)所述的频率、幅值和相位结合(2)所述的有关频谱校正的表达形式，对相位进行修正，获得相位的校正值；

获得步骤(1)所述的三个求和值，包括以下步骤：

(1.1)将待分析的离散时间信号设为x(t)，采样长度为2N，采样时间为2T；

(1.2)截取该信号[T／2，3T／2)区间上的部分作为x₁(t)，[0，T)区间上的部分作为x₂(t)，[T，2T)区间上的部分作为x₃(t)；

(1.3)则每段信号时长为T，长度为N，相邻信号的时间延迟为T／2；

(1.4)对每段时间信号分别求和，获得三个求和值：其中X_i表示编号为i的三段时间信号的求和值，i＝1，2，3；

获得步骤(2)所述的离散频谱低频成分的有关频谱校正的表达形式的方法如下：

(2.1).设三段信号的傅立叶变换为X_i(f)，i＝1，2，3，f表示频率；

(2.2).利用离散时间信号频谱的对称性，使用傅立叶变换公式，可获得离散频谱低频成分的频谱校正的表达形式为：

$\cos \left(\mathrm{\π}{\mathrm{Tf}}_0\right)=\frac{X_2\left(0\right)+X_3\left(0\right)}{{2X}_1\left(0\right)}$

$\tan \mathrm{\θ}=\frac{X_3\left(0\right)\cos \left(\mathrm{\π}{\mathrm{Tf}}_0\right)-X_1\left(0\right)}{X_3\left(0\right)\sin \left(\mathrm{\π}{\mathrm{Tf}}_0\right)}$

$A=\frac{X_3\left(0\right)}{W\left(f_0\right)\cos \left(\mathrm{\θ}\right)}$

式中，f₀表示离散时间信号中低频成分的频率校正值，X_i(0)分别表示三段信号频谱0Hz处的谱线值，其中i＝1，2，3，θ表示低频成分的相位，A表示低频成分的幅值校正值，W(f₀)表示窗函数频谱在频率f₀处的幅值；

(2.3).利用离散傅立叶变换公式的性质，可获得X_i(0)与X_i之间的关系为：X_i(0)=X_i，其中i＝1，2，3。

2.根据权利要求1所述的基于时间延迟的离散频谱低频成分的校正方法，其特征在于获得步骤(3)所述的离散频谱中低频成分的幅值、频率和相位的校正值的具体步骤如下：

根据X_i(0)和X_i相等，其中i＝1，2，3，可将X_i，i＝1，2，3结合频谱校正的表达形式，得到频率f₀、相位θ和幅值A，其中频率f₀和幅值A为所要求的校正值，相位θ经过修正公式后得到所要求的相位校正值。

3.根据权利要求2所述的基于时间延迟的离散频谱低频成分的校正方法，其特征在于所述的修正公式为θ′=θ+180f₀／f_s，其中，f_s为采样频率，θ′即为所要求的相位校正值。

4.根据权利要求1所述的基于时间延迟的离散频谱低频成分的校正方法，其特征在于适用于对称窗函数。