CN101997788A - 一种信号恢复的优化方法 - Google Patents

Abstract

本发明属于信号处理领域，尤其适用于应用在地球物理领域。一种信号恢复的优化方法，所述方法在对原始信号做双谱运算前，先对原始信号做预处理过程；所述预处理过程包括去线性背景、零均值化、周期拓展；完成预处理过程的信号再经过双谱运算过程后，实现了原始信号的完全恢复。本发明使基频f₀的估算既简便又准确，而且使基于双谱运算的信号恢复这一技术，能够适用于任意信号的信号恢复处理。本发明为信号恢复、图像降噪等处理提供了一个有用的工具。

Description

一种信号恢复的优化方法

技术领域

本发明属于信号处理领域，尤其适用于应用在地球物理、通讯工程、电子工程、军事等多个应用领域。

背景技术

现有技术中，基于双谱运算的信号重构，最重要的一个步骤是对基频率(或基波数)f₀的准确估算。现有的方法为：取双谱的对角线值B(f₁，f₂)，使|B(f，f)|不为0的最小f便是估算的f′₀。f′₀为f₀的初次估值，不一定精确，特别是当基频f₀与频率分辨率Δf之比不为整数时。估算的f′₀与真值f₀的误差范围为为减小估算值f′₀与真值f₀的误差，进一步搜索

区间内使|B(f，f)|达到最大的f，进行频率细化处理。对于实际数据，由于噪声的影响，|B(f，f)|不可能为0，采用第一个谱峰来代替不为0的|B(f，f)|。在进行频率细化处理的搜索计算时，先设定误差控制的精度ε。搜索计算采用迭代算法：

(1)令B＝|B(f′₀，f′₀)|，

计算|B(f，f)|；

(2)若B＜|B(f，f)|，则令B＝|B(f，f)|，f′₀＝f，若Δf＜ε，则结束，否则转回(1)；若B＞|B(f，f)|，则转入(3)。

(3)令B＝|B(f′₀，f′₀)|，计算|B(f，f)|；

若B＜|B(f，f)|，则令B＝|B(f，f)|，f′₀＝f，

若Δf＜ε，则结束，否则转回(1)；若B＞|B(f，f)|，则令

若Δf＜ε，则结束，否则转回(1)。

为减小估算值f′₀与真值f₀的误差，在

区间内，迭代搜索使|B(f，f)|达到最大的f，作为实际计算得到的基频f₀。

现有的技术中存在以下3方面的技术问题：

(1)通常情况下，信号的基频f₀与频率分辨率Δf之比不为整数，因此现有技术的方法估算的f₀，在信噪比较低时是个近似值(如表1)；即使基频f₀与频率分辨率Δf之比为整数，当信噪比较低时，上述方法估算的f₀也不能确保得到真值(如表2)。为使基频f₀估算的精度足够高，迭代运算的次数、计算量随之大大增加。

表1半基频f₀/2不为频率分辨率Δf的整数倍时不同信噪比信号估算的f′₀及其双谱值|B(f′₀，f′₀)|

表2半基频f₀/2为频率分辨率Δf的整数倍时不同信噪比信号估算的f′₀及其双谱值|B(f′₀，f′₀)|

(2)现有技术的方法，频率分辨率Δf越高，则基频f₀估算的精度越高，但对于实际信号，计算用的信号长度与信号的时间域采样率已经固定，也即频率分辨率Δf已经固定，因此，增加迭代运算的次数是提高估算精度的唯一途径。

(3)现有技术由于在方法上存在缺陷而带来实际应用上的局限。目前仅用于具有周期变化的信号的恢复处理，例如，仿真试验用的模型信号是具有周期变化的谐波信号，实际信号试验用的也是具有周期变化的缝纫机振动信号。

为了克服现有技术中基频f₀估算不准的缺陷，增强“应用双谱运算进行信号恢复”这一技术的适用性和实用性，申请人提出了本发明。

发明内容

本发明为了克服上述现有技术中存在的技术问题，研发了一种信号恢复的优化方法。本发明使基于双谱运算的信号重构技术适用于任意信号数据的重构，而且重构信号的保真度高。

本发明的核心内容是，通过信号处理方案、流程的改变，规避现有技术中方法的缺陷。具体地说，就是在对原始信号做双谱运算前，先对原始信号做3项预处理：去线性背景；零均值化；周期拓展。

其中，最关键之处是，对原始信号进行周期拓展处理。因为周期拓展，使信号的基频f₀与频率分辨率Δf相等。这样，不管原始信号的信噪比如何，基频f₀与频率分辨率Δf的比值始终是整数1，从而，确保了估算得到的基频f₀是真值。具体实现时，通过搜索双谱值|B(f，f)|的第一个谱峰时的频率f即可得到基频f₀的真值。为使后面应用式(14)推算谐波项系数a_i与相位角

时，方便、准确地搜索、应用双谱振幅谱B_y(f₀/2，f₀/2)和相位谱φ_y(f₀/2，f₀/2)的数值，根据经验将原始信号拓展成4个周期比较合适。

具体方法如下，

一种信号恢复的优化方法，所述方法在对地球物理勘探采集信号做双谱运算前，先对原始信号做预处理过程；所述预处理过程包括去线性背景、零均值化、周期拓展；完成预处理过程的信号再经过双谱运算过程后，实现了原始信号的完全恢复：完整描述一个信号的谐波信号基频、谐波项系数与相位等3项参数被准确求取。

所述方法包括如下步骤：

步骤1读取输入离散化后地球物理采集信号x(t)的数据：x(i)，i＝1，…，n；其中，n为信号x(t)的采样点数；

步骤2对信号x(t)进行预处理过程：

(1)去线性背景过程：

利用以下公式计算信号x(t)的线性背景斜率bg：

$\mathrm{bg}=\frac{x\left(n\right)-x\left(1\right)}{n-1}$

式中，x(n)为信号x(t)的最末一个数据，x(1)为信号x(t)的第一个数据，n为信号x(t)的采样点数；

然后利用以下公式消去信号的线性背景：

xb(i)＝x(i)-x(1)-bg×(i-1)

式中，xb(i)，i＝1，…，n为信号x(t)消去线性背景后的数据，x(i)，i＝1，…，n为原始信号x(t)的数据，x(1)为原始信号x(t)的第一个数据；

(2)零均值化过程：

对信号x(t)消去线性背景后的数据xb(i)，i＝1，…，n的各个样点数据xb(i)，求和并除以样点数n，得信号去线性背景后数据xb(i)，i＝1，…，n的平均值

$\stackrel{\_}{\mathrm{xb}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{xb}\left(i\right)/n$

式中，n为信号x(t)的采样点数；

去线性背景后的数据xb(i)，i＝1，…，n减去平均值

实现零均值化：

$\mathrm{xe}\left(i\right)=\mathrm{xb}\left(i\right)-\stackrel{\_}{\mathrm{xb}}$

式中，xe(i)，i＝1，…，n为信号x(t)去线性背景后又零均值化后的数据；

(3)周期拓展过程：

原信号是周期数为1的信号，将原信号首尾相连重复放置于后，重复放置K-1次，即得周期数为K的新信号；

周期拓展过程的C语言计算程序语句如下：

for(i＝0；i＜k；i++)for(j＝0；j＜n；j++)xn[j+i*k]＝xe[j]；

其中，xe[j]为原始信号x(t)去线性背景后又零均值化后数据xe(i)，i＝1，…，n的数组，xn[j+i*k]为周期拓展后信号xn(t)的数组，n为信号x(t)的采样点数，k为周期拓展后信号的周期数；

步骤3通过Fourier变换，计算周期拓展后信号xn(t)的频谱X(f)：

X(f)＝A(f)·e^jφ(f)

其中，A(f)和φ(f)分别为频谱X(f)的幅值和相位；其中，频谱X(f)的长度为k·n；其中，n为信号x(t)的采样点数，k为周期拓展后信号的周期数；

步骤4根据下式，通过三重相关计算信号xn(t)的双谱B(f₁，f₂)：

$B(f_1,f_2)=X\left(f_1\right)\·X\left(f_2\right)\·X^{*}(f_1+f_2)$

$=A\left(f_1\right)\·A\left(f_2\right)\·A(f_1+f_2)\·e^{j[\mathrm{\φ}\left(f_1\right)+\mathrm{\φ}\left(f_2\right)-\mathrm{\φ}(f_1+f_2)]}$

$=A(f_1,f_2)\·e^{\mathrm{j\φ}(f_1,f_2)}$

其中，X(f₁)、X(f₂)、X(f₁+f₂)分别是信号xn(t)在频率变量f₁，f₂，f₁+f₂上的频谱，X^*(f₁+f₂)是频谱X(f₁+f₂)的共轭，A(f₁，f₂)是双谱B(f₁，f₂)的幅值，φ(f₁，f₂)是双谱B(f₁，f₂)的相位；其中f₁，f₂是与时间变量t₁，t₂相对应的频率变量，t₁，t₂表示两个不同的时间变量；

这一过程的C语言计算程序语句如下：

for(i＝0；i＜km；i++){

  for(j＝0；j＜＝i；j++){

    if(i+j＞＝km)break；

    ba[i][j]＝a[i]*a[j]*a[i+j]；

    bf[i][j]＝f[i]+f[j]-f[i+j]；

  }

}

其中，数组a[i]、f[i]分别是信号xn(t)的频谱X(f)的振幅谱和相位谱数组，数组ba[i][j]、bf[i][j]分别是信号xn(t)的双谱B(f₁，f₂)的振幅谱A(f₁，f₂)和相位谱φ(f₁，f₂)数组；其中km信号xn(t)的频谱X(f)的长度的1/2，即km＝k·n/2+1；

步骤5利用对称性计算整个第一象限的双谱B(f₁，f₂)：

B(f₂，f₁)＝B(f₁，f₂)

其中，f₁，f₂的含义同上；

这一过程的C语言计算程序语句如下：

for(i＝0；i＜km；i++){

  for(j＝0；j＜＝i；j++){

    if(i+j＞＝km)break；

    ba[j][i]＝ba[i][j]；

    bf[j][i]＝bf[i][j]；

  }

}

其中，数组ba[i][j]、bf[i][j]分别是信号xn(t)的双谱B(f₁，f₂)的振幅谱数组和相位谱数组；其中km信号xn(t)的频谱X(f)的长度的1/2，即km＝k·n/2+1；

步骤6取信号xn(t)的双谱B(f₁，f₂)的模|B(f，f)|的第一个峰值时的f为估算的基频f₀；

步骤7在信号xn(t)上加一相位为0的半基频余弦信号，即y(t)＝xn(t)+cos(πf₀t)；

步骤8同步骤3，通过Fourier变换，计算信号y(t)的频谱Y(f)：

Y(f)＝A(f)·e^jφ(f)

其中，A(f)和φ(f)分别为频谱Y(f)的幅值和相位；频谱Y(f)的长度为k·n；其中，n为信号x(t)的采样点数，k为周期拓展后信号的周期数；

步骤9同步骤4，通过三重相关计算信号y(t)的双谱B_y(f₁，f₂)；

步骤10同步骤5，利用对称性计算整个第一象限的双谱B_y(f₁，f₂)；

步骤11设所要恢复的信号x(t)可以用谐波信号来描述，即

其中，f₀是谐波基频，已在步骤6求得；a_i，

i＝1，…，M，分别是谐波项系数、相位角，可以根据下式由信号y(t)的双谱B_y(f₁，f₂)的振幅谱A_y(f₁，f₂)与相位谱φ_y(f₁，f₂)求得：

$\left\{\begin{array}{c}{a}_1=8A_y(\frac{f_0}{2},\frac{f_0}{2})\\ a_i=\frac{8A_y[f_0,(i-1)f_0]}{a_1{a}_{i-1}}\end{array}\right.$ (i＝2，3，…，M)

(i＝2，3，…，M)

其中，M是谐波项数，可由下式求得：

M＝km/nf

其中，km为信号xn(t)的频谱X(f)的长度的1/2，即km＝k·n/2+1；nf是估算得到的基频f₀所对应的频率采样序号；

步骤12将基频f₀、谐波项系数a_i、谐波项相位角

谐波项数M，一起代入描述信号x(t)的谐波信号表达式：

得重构信号；

步骤13输出信号x(t)恢复结果。

本方法在步骤12后还包括进行线性背景补偿和均值补偿处理得恢复信号步骤；

所述均值补偿处理过程为：

对输出信号x(t)，加上步骤2中零均值化所求得的平均值

即 $x^{\′}\left(t\right)=x\left(t\right)+\stackrel{\_}{\mathrm{xb}};$

所述线性背景补偿过程为：

对均值补偿处理后的信号x′(t)，再加上步骤2中去线性背景时所求得的线性背景斜率bg：

x″(i)＝x′(i)+x(1)+bg×(i-1)

其中，x″(i)，i＝1，…，n为最终恢复得到的信号x″(t)的数据，x′(i)，i＝1，…，n为输出信号x(t)作了均值补偿处理后的数据，x(1)为原始信号的第一个数据，bg为原始信号的线性背景斜率。

本发明使基频f₀的估算既简便又准确，而且使基于双谱运算的信号恢复这一技术，能够适用于任意信号的信号恢复处理。本发明为信号恢复、图像降噪等处理提供了一个有用的工具。通过对现有技术文献中相同的谐波信号模型、地球物理勘探中的重力测量数据、地球物理勘探中的地震子波信号、地球物理勘探中的实际地震数据等多种信号数据，应用本发明研制的技术进行信号重构计算，基频f₀估算准确，重构信号的保真度高，方法技术的适用性和实用性强。

附图说明

图1为本发明方法流程框图；

图2为地球物理勘探中的重力测量信号的恢复处理结果；

图3为地球物理勘探中的地震子波信号恢复处理的结果；

图4为地球物理勘探中的原始地震信号剖面；

图5为应用本发明的方法恢复处理的地震信号剖面；

图6为放大100倍显示的附图4与附图5的差值。

下面将结合具体实施方式对各幅附图进行说明

具体实施方式

本发明所采用的方法流程详细说明：

第一步，读取离散化后地球物理信号x(t)的数据文件。

第二步，对信号x(t)进行去线性背景、零均值化、周期拓展等预处理。

第三步，通过Fourier变换，计算信号x(t)的频谱X(f)。

第四步，根据式(4)，通过三重相关计算信号x(t)的双谱B(f₁，f₂)。

第五步，利用对称性计算整个第一象限的双谱B(f₁，f₂)。

第六步，取|B(f，f)|的第一个峰值时的f为估算的基频f₀。

第七步，在原信号x(t)上加一相位为0的半基频余弦信号，即y(t)＝x(t)+cos(πf₀t)。

第八步，通过Fourier变换，计算信号y(t)的频谱Y(f)。

第九步，根据式(4)，通过三重相关计算信号y(t)的双谱B_y(f₁，f₂)。

第十步，利用对称性计算整个第一象限的双谱B_y(f₁，f₂)。

第十一步，根据式(14)计算谐波项系数a_i与相位角

第十二步，谐波项系数a_i与相位角代入式(5)，得重构信号。

第十三步，进行线性背景补偿和均值补偿处理得恢复信号。

第十四步，输出信号恢复结果。

在具体各个步骤中，依次为

(1)基于双谱运算信号恢复的方法原理

设x(t)是一均值为0的实信号，其Fourier变换为：

X(f)＝A(f)·e^jφ(f)                (1)

A(f)和φ(f)分别为频谱X(f)的幅值和相位。

x(t)的三阶累积量函数表达式为：

$C(t_1,t_2)={\∫}_{-\∞}^{+\∞}x\left(t\right)\·x(t+t_1)\·x(t+t_2)\mathrm{dt}---\left(2\right)$

t₁，t₂表示两个不同的时间变量。

C(t₁，t₂)的二维Fourier变换即为双谱B(f₁，f₂)

$B(f_1,f_2)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}C(t_1,t_2)\·e^{-j2\mathrm{\π}(f_1{t}_1+f_2{t}_2)}{\mathrm{dt}}_1{\mathrm{dt}}_2---\left(3\right)$

$=X\left(f_1\right)\·X\left(f_2\right)\·X^{*}(f_1+f_2)$

X^*(f₁+f₂)是频谱X(f₁+f₂)的共轭(注：下同)。其中f₁，f₂是与时间变量t₁，t₂相对应的频率变量。

进一步，可得x(t)的双谱表达式：

$B(f_1,f_2)=X\left(f_1\right)\·X\left(f_2\right)\·X^{*}(f_1+f_2)$

$=A\left(f_1\right)\·A\left(f_2\right)\·A(f_1+f_2)\·e^{j[\mathrm{\φ}\left(f_1\right)+\mathrm{\φ}\left(f_2\right)-\mathrm{\φ}(f_1+f_2)]}---\left(4\right)$

$=A(f_1,f_2)\·e^{\mathrm{j\φ}(f_1,f_2)}$

设所要恢复的信号x(t)可以用谐波信号来描述，即

根据Fourier变换的特性，由式(4)和(5)可以得到双谱振幅谱、相位谱与谐波项系数、相位角的关系式：

$A({\mathrm{if}}_0,{\mathrm{jf}}_0)=\frac{a_i{a}_j{a}_{i+j}}{8}$ (i＝1，2，…，M/2，j＝i，i+1，…，M-i)            (6a)

(i＝1，2，…，M/2，j＝i，i+1，…，M-i)            (6b)

方程式(6a)、(6b)中的未知数也即式(5)中的未知数，有基频f₀、谐波项系数a_i与相位角

(i＝2，3，…，M)。显然，(6a)、(6b)中未知数的个数大于方程数。

实际应用时须先对基频f₀进行准确地估算，这是本方法成功实现的关键，而且，

由式(5)可知，基频f₀估算的准确度，决定了恢复信号的保真度。

根据估算得到的基频f₀，在信号x(t)上加一相位为0的半基频余弦信号

z(t)＝cos(πf₀t)，即

y(t)＝x(t)+z(t)＝x(t)+cos(πf₀t)        (7)

由于y(t)的Fourier变换可以表示为：

Y(f)＝X(f)+Z(f)                         (8)

故y(t)的双谱可以表示为：

B_y(f₁，f₂)＝Y(f₁)·Y(f₂)·Y^*(f₁+f₂)

＝[X(f₁)+Z(f₁)]·[X(f₂)+Z(f₂)]·[X^*(f₁+f₂)+Z^*(f₁+f₂)]        (9)

对于式(9)，如果f₁＝f₂＝f₀/2，则根据Fourier变换的特性，有

从而有

$B_y(\frac{f_0}{2},\frac{f_0}{2})=[X\left(\frac{f_0}{2}\right)+Z\left(\frac{f_0}{2}\right)]\·[X\left(\frac{f_0}{2}\right)+Z\left(\frac{f_0}{2}\right)]\·[X^{*}\left(f_0\right)+X^{*}\left(f_0\right)]---\left(11a\right)$

如果f₁、f₂均不为f₀/2，且f₁+f₂≠f₀/2，则有

Z(f₁)＝Z(f₂)＝Z^*(f₁+f₂)＝0                                   (12)

从而，由式(12)、(9)可得

B_y(f₁，f₂)＝X(f₁)·X(f₂)·X^*(f₁+f₂)＝B_x(f₁，f₂)              (13)

即信号y(t)的双谱此时就是信号x(t)的双谱。从而可得所求谐波信号表达式x(t)中的谐波项系数a_i与相位角

求解的递推公式：

$\left\{\begin{array}{c}{a}_1=8A_y(\frac{f_0}{2},\frac{f_0}{2})\\ a_i=\frac{8A_y[f_0,(i-1)f_0]}{a_1{a}_{i-1}}\end{array}\right.$ (i＝2，3，…，M)                                                     (14a)

(i＝2，3，…，M)                                                     (14b)

由此，即实现了对信号x(t)的恢复。

根据上述技术流程，发明人在Windows操作系统下利用C++Builder，依照图1的程序框图，实现了采用本发明所述的方法做的任意信号的恢复处理，图中虚框部分为本发明的核心内容。

(1)地球物理勘探中的重力测量信号的恢复处理(附图2)。除了个别的单点异常外，原始信号几乎得到了完全的恢复。这些个别的单点，则反映了信号采集时的随机噪声的存在。说明本发明的信号恢复方法，可以实现对信号中随机噪声的消除。

(2)地球物理勘探中的地震子波信号的恢复处理(附图3)。地震子波的主频为202.4Hz。这一信号近乎脉冲信号，但本发明的方法也使它实现了信号的完全恢复，恢复信号与原始模型信号的均方误差为0.00076106。

(3)地球物理勘探中的实际地震记录信号的恢复处理(附图4～附图6)。从附图5与附图4的对比来看，原始地震记录信号得到了完全的恢复，但从附图6上来看，恢复信号与原始信号的差值是一些在“0”值附近来回震荡的随机噪声。说明应用本发明的方法可以使图像中的随机噪声有效地消除。

Claims (4)

1.一种信号恢复的优化方法，其特征在于，所述方法在对地球物理采集信号做双谱运算前，先对原始信号做预处理过程；所述预处理过程包括去线性背景、零均值化、周期拓展；完成预处理过程的信号再经过双谱运算过程后，实现了原始信号的完全恢复，实现完整描述一个信号的谐波信号基频、谐波项系数与相位参数。

2.根据权利要求1所述的一种信号恢复的优化方法，其特征在于，所述方法包括如下步骤：

步骤1读取输入离散化后地球物理采集信号x(t)的数据：x(i)，i＝1，…，n；其中，n为信号x(t)的采样点数；

步骤2对信号x(t)进行预处理过程：

(1)去线性背景过程：

利用以下公式计算信号x(t)的线性背景斜率bg：

$\mathrm{bg}=\frac{x\left(n\right)-x\left(1\right)}{n-1}$

式中，x(n)为信号x(t)的最末一个数据，x(1)为信号x(t)的第一个数据，n为信号x(t)的采样点数；

然后利用以下公式消去信号的线性背景：

xb(i)＝x(i)-x(1)-bg×(i-1)

式中，xb(i)，i＝1，…，n为信号x(t)消去线性背景后的数据，x(i)，i＝1，…，n为原始信号x(t)的数据，x(1)为原始信号x(t)的第一个数据；

(2)零均值化过程：

对信号x(t)消去线性背景后的数据xb(i)，i＝1，…，n的各个样点数据xb(i)，求和并除以样点数n，得信号去线性背景后数据xb(i)，i＝1，…，n的平均值

$\stackrel{\_}{\mathrm{xb}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{xb}\left(i\right)/n$

式中，n为信号x(t)的采样点数；

去线性背景后的数据xb(i)，i＝1，…，n减去平均值

实现零均值化：

$\mathrm{xe}\left(i\right)=\mathrm{xb}\left(i\right)-\stackrel{\_}{\mathrm{xb}}$

式中，xe(i)，i＝1，…，n为信号x(t)去线性背景后又零均值化后的数据；

(3)周期拓展过程：

原信号是周期数为1的信号，将原信号首尾相连重复放置于后，重复放置K-1次，即得周期数为K的新信号；

步骤3通过Fourier变换，计算周期拓展后信号xn(t)的频谱X(f)：

X(f)＝A(f)·e^jφ(f)

其中，A(f)和φ(f)分别为频谱X(f)的幅值和相位；其中，频谱X(f)的长度为k·n；其中，n为信号x(t)的采样点数，k为周期拓展后信号的周期数；

步骤4通过三重相关计算得到信号xn(t)的双谱B(f₁，f₂)：

$B(f_1,f_2)=X\left(f_1\right)\·X\left(f_2\right)\·X^{*}(f_1+f_2)$

$=A\left(f_1\right)\·A\left(f_2\right)\·A(f_1+f_2)\·e^{j[\mathrm{\φ}\left(f_1\right)+\mathrm{\φ}\left(f_2\right)-\mathrm{\φ}(f_1+f_2)]}$

$=A(f_1,f_2)\·e^{\mathrm{j\φ}(f_1,f_2)}$

其中，X(f₁)、X(f₂)、X(f₁+f₂)分别是信号xn(t)在频率变量f₁，f₂，f₁+f₂上的频谱，X^*(f₁+f₂)是频谱X(f₁+f₂)的共轭，A(f₁，f₂)是双谱B(f₁，f₂)的幅值，φ(f₁，f₂)是双谱B(f₁，f₂)的相位；其中f₁，f₂是与时间变量t₁，t₂相对应的频率变量，t₁，t₂表示两个不同的时间变量；

步骤5利用对称性计算整个第一象限的双谱B(f₁，f₂)：

B(f₂，f₁)＝B(f₁，f₂)

其中，f₁，f₂的含义同上；

步骤6取信号xn(t)的双谱B(f₁，f₂)的模|B(f，f)|的第一个峰值时的f为估算的基频f₀；

步骤7在信号xn(t)上加一相位为0的半基频余弦信号，即y(t)＝xn(t)+cos(πf₀t)；

步骤8同步骤3，通过Fourier变换，计算信号y(t)的频谱Y(f)：

Y(f)＝A(f)·e^jφ(f)

其中，A(f)和φ(f)分别为频谱Y(f)的幅值和相位；频谱Y(f)的长度为k·n；其中，n为信号x(t)的采样点数，k为周期拓展后信号的周期数；

步骤9同步骤4，通过三重相关计算信号y(t)的双谱B_y(f₁，f₂)；

步骤10同步骤5，利用对称性计算整个第一象限的双谱B_y(f₁，f₂)；

步骤11设所要恢复的信号x(t)可以用谐波信号来描述，即

其中，f₀是谐波基频，已在步骤6求得；a_i，

i＝1，…，M，分别是谐波项系数、相位角，可以根据下式由信号y(t)的双谱B_y(f₁，f₂)的振幅谱A_y(f₁，f₂)与相位谱φ_y(f₁，f₂)求得：

$\left\{\begin{array}{c}{a}_1=8A_y(\frac{f_0}{2},\frac{f_0}{2})\\ a_i=\frac{8A_y[f_0,(i-1)f_0]}{a_1{a}_{i-1}}\end{array}\right.$ (i＝2，3，…，M)

(i＝2，3，…，M)

其中，M是谐波项数，可由下式求得：

M＝km/nf

其中，km为信号xn(t)的频谱X(f)的长度的1/2，即km＝k·n/2+1；nf是估算得到的基频f₀所对应的频率采样序号；

步骤12将基频f₀、谐波项系数a_i、谐波项相位角

谐波项数M，一起代入描述信号x(t)的谐波信号表达式：

得重构信号；

步骤13输出信号x(t)恢复结果。

3.根据权利要求2所述的一种信号恢复的优化方法，其特征在于，所述方法包括如下步骤：

本方法在步骤12后还包括进行线性背景补偿和均值补偿处理得恢复信号步骤；

所述均值补偿处理过程为：

对输出信号x(t)，加上步骤2中零均值化所求得的平均值

即 $x^{\′}\left(t\right)=x\left(t\right)+\stackrel{\_}{\mathrm{xb}};$

所述线性背景补偿过程为：

对均值补偿处理后的信号x′(t)，再加上步骤2中去线性背景时所求得的线性背景斜率bg：

x″(i)＝x′(i)+x(1)+bg×(i-1)

其中，x″(i)，i＝1，…，n为最终恢复得到的信号x″(t)的数据，x′(i)，i＝1，…，n为输出信号x(t)作了均值补偿处理后的数据，x(1)为原始信号的第一个数据，bg为原始信号的线性背景斜率。

4.根据权利要求2所述的一种信号恢复的优化方法，其特征在于，所述方法的步骤2中周期拓展过程：采用将原始信号数据拓展成4个周期。

Images