CN102570984B - 一种用于电学层析成像***的多频递推解调方法 - Google Patents

Abstract

本发明涉及一种用于电学层析成像***的多频递推解调方法，其特征在于：其包括下列步骤：(1)根据激励信号的m个分量频率和***采样频率建立初始方程和递推方程；(2)以激励信号零相位时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2m个采样点数据按照步骤(1)中的初始方程计算出数据矩阵B和状态矩阵P的初始值；(3)将其初始值及新增采样点数据代入步骤(1)中的递推方程，以更新B和P；(4)判断解调结果是否满足***精度要求，若不满足则返回步骤(3)，若满足则结束解调，并输出测量信号的m个频率分量的幅值和相位数据。上述递推解调方法最少只需2m个采样点数据，且增加采样点数可提高解调精度和抗噪声能力。

Description

一种用于电学层析成像***的多频递推解调方法

技术领域

本发明涉及一种用于电学层析成像***的多频递推解调方法，属于分布参数测量领域。

背景技术

电学层析成像（Electrical Tomography）通过对被测对象施加电激励，并检测其边界值的变化，利用特定数学手段反演被测对象内部的电学参数分布，从而得到对象内部的物质分布。与其它层析成像技术相比，电学层析成像具有无辐射、响应速度快、价格低廉等优势。

由于电学层析成像多应用于工业多相流检测和医学检测等领域，被测对象通常表现为相对复杂电学特性，且其随激励频率的变化而变化。因此，为获得被测对象内部电学参数更加丰富的频谱特性，越来越多的电学成像***采用多种频率正弦信号相叠加的混频激励方式。然而，混频激励的方式也给测量带来了更多的困难，尤其体现在解调方式上。

为实现多频率信号解调的相敏解调，通常将测量信号的每一个频率分量与激励信号的对应分量分别解调以得到测量信号在不同频率分量下幅值和相位信息。传统的模拟多频相敏解调的方式需要与激励频率分量数量相同的相敏解调器，大大的增加了硬件电路规模和复杂度，降低了***的实时性。随着数字技术的高速发展，基于数字处理器的数字相敏解调方法以其良好的测量实时性和简洁的硬件电路结构而受到研究者越来越多的关注。

数字相敏解调方法首先利用高速模数转换器对被测信号进行采样，然后利用高性能数字信号处理器件，如FPGA、DSP等，采用数值计算的方法提取被测信号的幅值和相位信息。目前最为常用的数字相敏解调方法为正交序列解调方法，应用到多频电学层析成像***中时，要求采样序列长度必须覆盖混频激 励信号的整数个周期（由最低频率分量的信号周期决定）才能解调出测量信号在所有频率分量下的幅值和相位信息，这在很大程度上降低了解调方法的灵活性，同时也限制了解调速度的进一步提高。

发明内容

本发明的目的在于提供一种用于电学层析成像***的多频递推解调方法，可在少于一个完整混频信号周期的时间内得到精度较高的解调结果，且随着代入递推过程的采样点数的增加，可提高解调结果的抗噪声性能。

本发明所提供的一种用于电学层析成像***的多频递推解调方法，包括下列步骤：

步骤一、根据***激励信号的m个分量频率f₁，f₂，…，f_m以及采样频率f_s建立初始方程和递推方程：

初始方程： $\left\{\begin{array}{c}P(2m-1)={[V{(2m-1)}^H\·V(2m-1)]}^{-1}\\ B_{2m-1}=P{(2m-1)}^{-1}\·V{(2m-1)}^H\·X(2m-1)\end{array}\right.---\left(1\right)$

递推方程： $\left\{\begin{array}{c}P(k+1)=P\left(k\right)-\frac{P\left(k\right){V_{k+1}}^H{V}_{k+1}P\left(k\right)}{1+V_{k+1}P\left(k\right){V_{k+1}}^H}\\ B_{k+1}=B_k+P(k+1){V_{k+1}}^k(x_{k+1}-V_{k+1}{B}_k)\end{array}\right.,k\≥2m-1---\left(2\right)$

在方程(1)和方程(2)中，k为采样点序数且k≥0；x_k为第k个采样点数据；P(k)为状态矩阵；B_k为数据矩阵；V_k为由激励信号的各个分量频率和***采样频率唯一确定的常数向量，其具体形式为：

$V_k=\left[\begin{array}{cccccccc}{e}^{-k\·j2\mathrm{\π}{f}_1/f_s}& e^{-k\·j2\mathrm{\π}{f}_2/f_s}& ...& e^{-k\·j2\mathrm{\π}{f}_m/f_s}& e^{k\·j2\mathrm{\π}{f}_m/f_s}& ...& e^{k\·j2\mathrm{\π}{f}_2/f_s}& e^{k\·j2\mathrm{\π}{f}_1/f_s}\end{array}\right]---\left(3\right)$

V(2m-1)为由前2m个V_k构成的辅助矩阵，

V(2m-1)＝[V₀ V₁ … V_2m-1]^T (4)

X(2m-1)为由前2m个采样点数据构成的测量向量，

X(2m-1)＝[x₀ x₁ … x_2m-1]^T (5)

步骤二、以含有多个频率分量的激励信号零相位时刻为起始时刻,根据测量信号在紧邻起始时刻之后的2m个采样点数据，按照步骤一中的初始方 程(1)计算出包含测量信号各频率分量的幅值与相位信息的数据矩阵B和状态矩阵P的初始值B_2m-1和P(2m-1)；

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2m)个采样点数据x_k代入到步骤一中建立的递推方程(2)中，逐步更新数据矩阵B和状态矩阵P的值；

步骤四、判断解调结果是否满足***精度要求，若不满足，则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果计算出测量信号各频率分量的幅值和相位

$\left\{\begin{array}{c}\widehat{A_1}=2\left|\widehat{b_1}\right|\\ \widehat{{\mathrm{\θ}}_1}=\arctan [\mathrm{Im}\left(\widehat{b_1}\right)/\mathrm{Re}\left(\widehat{b_1}\right)]\\ \widehat{A_2}=2\left|\widehat{b_2}\right|\\ \widehat{{\mathrm{\θ}}_2}=\arctan [\mathrm{Im}\left(\widehat{b_2}\right)/\mathrm{Re}\left(\widehat{b_2}\right)]\\ .\\ .\\ .\\ \widehat{A_m}=2\left|\widehat{b_m}\right|\\ \widehat{{\mathrm{\θ}}_m}=\arctan [\mathrm{Im}\left(\widehat{b_m}\right)/\mathrm{Re}\left(\widehat{b_m}\right)]\end{array}\right..---\left(6\right)$

本发明与现有技术相比的优点在于：可先利用从紧邻起始时刻起的2m个（m为混频激励信号的频率分量个数）测量信号采样点数据得到初步的解调结果，之后根据***的测量精度和实时性等要求适当增加采样点数，通过递推的方法来提高解调精度和抗噪声能力；在无噪声的情况下，该方法不需要采用完整混频信号周期的采样点即可得到准确的解调结果，在存在噪声的情况下，亦可获得较为理想的解调结果，具有很好的实时性和灵活性。

附图说明

图1为本发明提供的递推解调方法实施过程的流程图；

图2(a)为仿真实验中测量信号无噪声时100kHz频率分量幅值解调结果；

图2(b)为仿真实验中测量信号无噪声时100kHz频率分量相位解调结果；

图2(c)为仿真实验中测量信号无噪声时200kHz频率分量幅值解调结果；

图2(d)为仿真实验中测量信号无噪声时200kHz频率分量相位解调结果；

图2(e)为仿真实验中测量信号无噪声时300kHz频率分量幅值解调结果；

图2(f)为仿真实验中测量信号无噪声时300kHz频率分量相位解调结果；

图2(g)为仿真实验中测量信号含噪声时100kHz频率分量幅值解调结果；

图2(h)为仿真实验中测量信号含噪声时100kHz频率分量相位解调结果；

图2(i)为仿真实验中测量信号含噪声时200kHz频率分量幅值解调结果；

图2(j)为仿真实验中测量信号含噪声时200kHz频率分量相位解调结果；

图2(k)为仿真实验中测量信号含噪声时300kHz频率分量幅值解调结果；

图2(l)为仿真实验中测量信号含噪声时300kHz频率分量相位解调结果。

具体实施方式

本发明，即一种用于电学层析成像***的多频递推解调方法，包括下列步骤：

步骤一、根据***激励信号的m个分量频率f₁，f₂，…，f_m以及采样频率f_s建立多初始方程和递推方程。

假设多频激励电学层析成像***激励信号的分量频率分别为f₁，f₂，…，f_m，测量信号中对应频率分量的幅值分别为A₁，A₂，…，A_m，若采样频率f_s，则测量信号可表示为：

$x_k=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{A}_i\cos (2\mathrm{\πk}{f}_i/f_s+{\mathrm{\θ}}_i)---\left(1\right)$

其中，k为采样点序数，θ_i为测量信号的第i个频率分量的相位。根据欧拉公式，x_k可表示为：

$x_k=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(\frac{A_i}{2}{e}^{-j{\mathrm{\θ}}_i}\·e^{-2\mathrm{\πk}{f}_i/f_s}+\frac{A_i}{2}{e}^{j{\mathrm{\θ}}_i}\·e^{2\mathrm{\πk}{f}_i/f_s})---\left(2\right)$

写成矩阵的形式为：

$x_k=\left[\begin{array}{cccccccc}{e}^{-2\mathrm{\πk}{f}_1/f_s}& e^{-2\mathrm{\πk}{f}_2/f_s}& ...& e^{-2\mathrm{\πk}{f}_m/f_s}& e^{2\mathrm{\πk}{f}_m/f_s}& ...& e^{-2\mathrm{\πk}{f}_2/f_s}& e^{2\mathrm{\πk}{f}_1/f_s}\end{array}\right]\left[\begin{array}{c}\frac{A_1}{2}{e}^{-j{\mathrm{\θ}}_1}\\ \frac{A_2}{2}{e}^{-j{\mathrm{\θ}}_2}\\ .\\ .\\ .\\ \frac{A_k}{2}{e}^{-j{\mathrm{\θ}}_m}\\ \frac{A_k}{2}{e}^{j{\mathrm{\θ}}_m}\\ .\\ .\\ .\\ \frac{A_2}{2}{e}^{j{\mathrm{\θ}}_2}\\ \frac{A_1}{2}{e}^{j{\mathrm{\θ}}_1}\end{array}\right]---\left(3\right)$

因此，考虑构造解调辅助矩阵V和含有测量信号幅值与相位信息的数据矩阵B

$B={\left[\begin{array}{cccccccc}\frac{A_1}{2}{e}^{-j{\mathrm{\θ}}_1}& \frac{A_2}{2}{e}^{-j{\mathrm{\θ}}_2}& ...& \frac{A_m}{2}{e}^{-j{\mathrm{\θ}}_m}& \frac{A_m}{2}{e}^{j{\mathrm{\θ}}_m}& ...& \frac{A_2}{2}{e}^{j{\mathrm{\θ}}_2}& \frac{A_1}{2}{e}^{j{\mathrm{\θ}}_1}\end{array}\right]}^T---\left(5\right)$

则测量信号的采样序列与解调辅助矩阵之间的关系可以表示为：

X＝[x₀ x₁ … x_k …]^T＝VB (6)

其中

$\begin{array}{c}{V}_0=\left[\begin{array}{cccccccc}{e}^{-0\·j2\mathrm{\π}{f}_1/f_s}& e^{-0\·j2\mathrm{\π}{f}_2/f_s}& ...& e^{-0\·j2\mathrm{\π}{f}_m/f_s}& e^{0\·j2\mathrm{\π}{f}_m/f_s}& ...& e^{0\·j2\mathrm{\π}{f}_2/f_s}& e^{0\·j2\mathrm{\π}{f}_1/f_s}\end{array}\right]\\ V_1=\left[\begin{array}{cccccccc}{e}^{-1\·j2\mathrm{\π}{f}_1/f_s}& e^{-1\·j2\mathrm{\π}{f}_2/f_s}& ...& e^{-1\·j2\mathrm{\π}{f}_m/f_s}& e^{1\·j2\mathrm{\π}{f}_m/f_s}& ...& e^{1\·j2\mathrm{\π}{f}_2/f_s}& e^{1\·j2\mathrm{\π}{f}_1/f_s}\end{array}\right]\\ .\\ .\\ .\\ V_k=\left[\begin{array}{cccccccc}{e}^{-k\·j2\mathrm{\π}{f}_1/f_s}& e^{-k\·j2\mathrm{\π}{f}_2/f_s}& ...& e^{-k\·j2\mathrm{\π}{f}_m/f_s}& e^{k\·j2\mathrm{\π}{f}_m/f_s}& ...& e^{k\·j2\mathrm{\π}{f}_2/f_s}& e^{k\·j2\mathrm{\π}{f}_1/f_s}\end{array}\right]\\ .\\ .\\ .\end{array}---\left(7\right)$

$V\left(0\right)=V_0,V\left(1\right)=\left[\begin{array}{c}{V}_0\\ V_1\end{array}\right]=\left[\begin{array}{c}V\left(0\right)\\ V_1\end{array}\right],...,V\left(k\right)=\left[\begin{array}{c}V(k-1)\\ V_k\end{array}\right],...---\left(8\right)$

X(0)＝[x₀],X₁＝[x₀ x₁]^T＝[X(0) x₁]^T,…,X_k＝[X(k-1) x_k]^T,… (9)

那么方程(6)可以表示成另一种形式，即

X(k)＝V(k)B (10)

如果V(k)是方阵且其逆存在，则可通过

B＝V(k)^-1X(k) (11)

直接计算出数据矩阵B，进而解调出测量信号的幅值和相位。然而，V(k)是一个k+1行2m列的矩阵，除k+1＝2m的情况下均不能用方程(11)求解B。

考虑到矩阵V(k)的特点，重新构造求解数据矩阵B的公式，当k≥2m-1时其可表示为：

B_k＝[V(k)^HV(k)]^-1V(k)^HX(k) (12)

令

P(k)＝[V(k)^HV(k)]^-1 (13)

则

$\begin{array}{c}P(k+1)={\left[V{(k+1)}^{H}V(k+1)\right]}^{-1}\\ ={\left({\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]}^H\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]\right)}^{-1}\\ ={\left(\left[\begin{array}{cc}V{\left(k\right)}^H& {V_{k+1}}^H\end{array}\right]\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]\right)}^{-1}\\ ={(V{\left(k\right)}^{H}V\left(k\right)+{V_{k+1}}^H{V}_{k+1})}^{-1}\\ ={(P{\left(k\right)}^{-1}+V_{k+1}^H{V}_{k+1})}^{-1}\end{array}---\left(14\right)$

根据矩阵求逆的定理，有

(A+BCD)^-1＝A^-1-A^-1B(C^-1+DA^-1B)^-1DA^-1 (15)

令A＝P(k)^-1，C＝1，D＝V_k+1，可得

P(k+1)＝(A+B·1·B^H)^-1

＝A^-1-A^-1B(1+B^HA^-1B)^-1B^HA^-1 (16)

＝P(k)-P(k)V_k+1 ^H(1+V_k+1P(k)V_k+1 ^H)^-1V_k+1P(k)

因此在方程(12)中，若用第k次计算的结果更新第k+1次计算的相关参数，则

$\begin{array}{c}{B}_{k+1}=P(k+1)V{(k+1)}^{H}X(k+1)\\ =P(k+1){\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]}^H\left[\begin{array}{c}X\left(k\right)\\ x_{k+1}\end{array}\right]\\ =P(k+1)\left[\begin{array}{cc}{V\left(k\right)}^H& {V_{k+1}}^H\end{array}\right]\left[\begin{array}{c}X\left(k\right)\\ x_{k+1}\end{array}\right]\\ =P(k+1)[{V\left(k\right)}^{H}X\left(k\right)+{V_{k+1}}^H{x}_{k+1}]\\ =P(k+1)P{\left(k\right)}^{-1}{B}_k+P(k+1){V_{k+1}}^H{x}_{k+1}\end{array}---\left(17\right)$

将式(14)变换为

P(k+1)^-1＝P(k)^-1+V_k+1 ^HV_k+1 (18)

从而，

P(k)^-1＝P(k+1)^-1-V_k+1 ^HV_k+1 (19)

将式(19)代入式(17)，可得

B_k+1＝P(k+1)[P(k+1)^-1-V_k+1 ^HV_k+1]B_k+P_k+1V_k+1 ^Hx_k+1

＝B_k-P(k+1)V_k+1 ^HV_k+1B_k+P(k+1)V_k+1 ^Hx_k+1 (20)

＝B_k+P(k+1)V_k+1 ^H(x_k+1-V_k+1B_k)

由式(16)和(20)，可知递推方程为

$\left\{\begin{array}{c}P(k+1)=P\left(k\right)-\frac{P\left(k\right){V_{k+1}}^H{V}_{k+1}P\left(k\right)}{1+V_{k+1}P\left(k\right){V_{k+1}}^H}\\ B_{k+1}=B_k+P(k+1){V_{k+1}}^k(x_{k+1}-V_{k+1}{B}_k)\end{array}\right.,k\≥2m-1---\left(21\right)$

而其初始值可由初始方程

$\left\{\begin{array}{c}P(2m-1)={[V{(2m-1)}^H\·V(2m-1)]}^{-1}\\ B_{2m-1}=P{(2m-1)}^{-1}\·V{(2m-1)}^H\·X(2m-1)\end{array}\right.---\left(22\right)$

计算得到。

步骤二、以含有多个频率分量的激励信号零相位时刻为起始时刻,根据测量信号在紧邻起始时刻之后的2m个采样点数据，按照步骤一中的初始方程(22)计算出包含测量信号各频率分量的幅值与相位信息的数据矩阵B和状态矩阵P的初始值B_2m-1和P(2m-1)；

令紧邻起始时刻之后的2m测量信号采样点的值为x₀，x₁，…，x_2m-1，根据方程(12)和(13)计算得到：

$\left\{\begin{array}{c}P(2m-1)={\left[V{(2m-1)}^{H}V(2m-1)\right]}^{-1}\\ B_{2m-1}=P(2m-1)V{(2m-1)}^{H}X(2m-1)={\left[\begin{array}{cccc}{b}_{11}& b_{12}& ...& b_{1\left(2m\right)}\end{array}\right]}^T\end{array}\right.---\left(23\right)$

此时，可根据初始数据矩阵B₁计算出测量信号的幅值和相角：

$\left\{\begin{array}{c}\widehat{A_{11}}=2\left|\widehat{b_{11}}\right|\\ \widehat{{\mathrm{\θ}}_{11}}=\arctan [\mathrm{Im}\left(\widehat{b_{11}}\right)/\mathrm{Re}\left(\widehat{b_{11}}\right)]\\ \widehat{A_{12}}=2\left|\widehat{b_{12}}\right|\\ \widehat{{\mathrm{\θ}}_{12}}=\arctan [\mathrm{Im}\left(\widehat{b_{12}}\right)/\mathrm{Re}\left(\widehat{b_{12}}\right)]\\ .\\ .\\ .\\ \widehat{A_{1m}}=2\left|\widehat{b_{1m}}\right|\\ \widehat{{\mathrm{\θ}}_{1m}}=\arctan [\mathrm{Im}\left(\widehat{b_{1m}}\right)/\mathrm{Re}\left(\widehat{b_{1m}}\right)]\end{array}\right..---\left(24\right)$

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2m)个采样点数据x_k代入到步骤一中建立的递推方程(21)中，逐步更新数据矩阵B和状态矩阵P的值；

步骤四、判断解调结果是否满足***精度要求，若不满足，则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果计算出测量信号各频率分量的幅值和相位

$\left\{\begin{array}{c}\widehat{A_1}=2\left|\widehat{b_1}\right|\\ \widehat{{\mathrm{\θ}}_1}=\arctan [\mathrm{Im}\left(\widehat{b_1}\right)/\mathrm{Re}\left(\widehat{b_1}\right)]\\ \widehat{A_2}=2\left|\widehat{b_2}\right|\\ \widehat{{\mathrm{\θ}}_2}=\arctan [\mathrm{Im}\left(\widehat{b_2}\right)/\mathrm{Re}\left(\widehat{b_2}\right)]\\ .\\ .\\ .\\ \widehat{A_m}=2\left|\widehat{b_m}\right|=2\left|\widehat{b_m}\right|\\ \widehat{{\mathrm{\θ}}_m}=\arctan [\mathrm{Im}\left(\widehat{b_m}\right)/\mathrm{Re}\left(\widehat{b_m}\right)]\end{array}\right..---\left(25\right)$

下面结合附图和实施例对本发明做进一步详细说明。

在Matlab计算软件中对本发明所提供的一种用于电学层析成像***的递推解调方法进行仿真实验。实验条件如下：

（1）测量信号为三个理想正弦信号分量的叠加，正弦分量一的幅值为1、相位为60°、频率为100kHz；正弦分量二的幅值为5、相位为45°、 频率为200kHz；正弦分量三的幅值为2、相位为30°、频率为300kHz；***的采样率为3MHz，每个混合信号周期内有30个采样点；

（2）在实验条件（1）中的理想混频测量信号中加入均值为0、标准差为0.01的随机噪声。

实验条件（1）下测量信号中三个频率分量幅值和相位的递推解调结果分别如图2中(a)-(f)所示，实验条件（2）下测量信号中三个频率分量的幅值和相位的递推解调结果分别如图2中(g)、(l)所示。

由仿真实验结果可以看出，在测量信号为3个理想正弦混频信号的情况下，本发明所提供的递推解调方法利用初始的6个采样点即可解调出各个频率分量幅值和相位的准确值，且随着采样点数的增多解调结果依然准确；而在测量信号含有3个频率分量的混频信号且含有随机噪声的情况下，本发明所提供的递推解调方法利用初始的6个采样点即可解调出精度较高的各频率分量的幅值和相位值，其相对误差均小于1%，且随着采样点数的增多解调相对误差呈现震荡衰减的趋势。

仿真实验验证了本发明所提供多频递推解调方法的良好效果。

以上对本发明及其实施方式的描述，并不局限于此，附图中所示仅是本发明的实施方式之一。在不脱离本发明创造宗旨的情况下，不经创造性地设计出与该技术方案类似的结构或实施例，均属本发明保护范围。

Claims (1)

1.一种用于电学层析成像***的多频递推解调方法，其特征在于：其包括下列步骤：

步骤一、根据***激励信号的m个分量频率f₁，f₂，…，f_m以及采样频率f_s建立初始方程和递推方程：

初始方程： $\left\{\begin{array}{c}P(2m-1)={[V{(2m-1)}^H\·V(2m-1)]}^{-1}\\ B_{2m-1}=P{(2m-1)}^{-1}\·V{(2m-1)}^H\·X(2m-1)\end{array}\right.---\left(1\right)$

递推方程： $\left\{\begin{array}{c}P(k+1)=P\left(k\right)-\frac{P\left(k\right){V_{k+1}}^H{V}_{k+1}P\left(k\right)}{1+V_{k+1}P\left(k\right){V_{k+1}}^H}\\ B_{k+1}=B_k+P(k+1){V_{k+1}}^k(x_{k+1}-V_{k+1}{B}_k)\end{array}\right.,k\≥2m-1---\left(2\right)$

在方程(1)和方程(2)中，k为采样点序数且k≥0；x_k为第k个采样点数据；P(k)为状态矩阵；B_k为数据矩阵；V_k为由激励信号的各个分量频率和***采样频率唯一确定的常数向量，其具体形式为：

$V_k=\left[\begin{array}{cccccccc}{e}^{-k\·j2\mathrm{\π}{f}_1/f_s}& e^{-k\·j2\mathrm{\π}{f}_2/f_s}& ...& e^{-k\·j2\mathrm{\π}{f}_m/f_s}& e^{k\·j2\mathrm{\π}{f}_m/f_s}& ...& e^{k\·j2\mathrm{\π}{f}_2/f_s}& e^{k\·j2\mathrm{\π}{f}_1/f_s}\end{array}\right]---\left(3\right)$

V(2m-1)为由前2m个V_k构成的辅助矩阵，

V(2m-1)＝[V₀ V₁ … V_2m-1]^T (4)

X(2m-1)为由前2m个采样点数据构成的测量向量，

X(2m-1)＝[x₀ x₁ … x_2m-1]^T (5)

步骤二、以含有多个频率分量的激励信号零相位时刻为起始时刻,根据测量信号在紧邻起始时刻之后的2m个采样点数据，按照步骤一中的初始方程(1)计算出包含测量信号各频率分量的幅值与相位信息的数据矩阵B和状态矩阵P的初始值B_2m-1和P(2m-1)；

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2m)个采样点数据x_k代入到步骤一中建立的递推方程(2)中，逐步更新数据矩阵B和状态矩阵P的值；

步骤四、判断解调结果是否满足***精度要求，若不满足，则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果计算出测量信号各频率分量的幅值和相位

$\left\{\begin{array}{c}\widehat{A_1}=2\left|\widehat{b_1}\right|\\ \widehat{{\mathrm{\θ}}_1}=\arctan [\mathrm{Im}\left(\widehat{b_1}\right)/\mathrm{Re}\left(\widehat{b_1}\right)]\\ \widehat{A_2}=2\left|\widehat{b_2}\right|\\ \widehat{{\mathrm{\θ}}_2}=\arctan [\mathrm{Im}\left(\widehat{b_2}\right)/\mathrm{Re}\left(\widehat{b_2}\right)]\\ .\\ .\\ .\\ \widehat{A_m}=2\left|\widehat{b_m}\right|\\ \widehat{{\mathrm{\θ}}_m}=\arctan [\mathrm{Im}\left(\widehat{b_m}\right)/\mathrm{Re}\left(\widehat{b_m}\right)]\end{array}\right..---\left(6\right)$