一种残差迭代的Prony低频振荡分析方法

摘要

本发明公开了一种残差迭代的Prony低频振荡分析方法，采用迭代的方法来减小信号模型的阶数，提高解算特征方程系数求解方程组的数值稳定性。具体包括以下步骤：读入滤波后的测量数据；用奇异值分解法计算信号模型实际阶数；用残差迭代法计算特征方程系数；计算频率和衰减因子；计算幅值和初相。发明通过残差迭代的Prony法分析低频振荡，可以减小扩展Prony分析时的信号模型的阶数，提高解算特征方程系数求解方程组的数值稳定性。在所取的信号模型阶数m相同的条件下，本发明比传统Prony法计算出的低频振荡信息的结果精度高。

主权项

一种残差迭代的Prony低频振荡分析方法，其特征在于：包括以下步骤：A、读入滤波后的测量数据读入测量数据，并对测量数据滤波去掉高频信号及噪声；B、用奇异值分解法计算信号模型实际阶数根据测量数据形成扩展阶矩阵X_e，所取信号模型阶数p_e大于信号模型的实际阶数p，p_e取值为满足[N/4]<p_e<[N/2]的整数，然后采用奇异值分解法计算扩展阶矩阵X_e的有效秩r，即得到信号模型的实际阶数p＝r；扩展阶矩阵X_e为：$X_e=\left[\begin{array}{cccc}x(p_e-1)& x(p_e-2)& ...& x\left(0\right)\\ x\left(p_e\right)& x(p_e-1)& ...& x\left(1\right)\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ x(N-2)& x(N-3)& ...& x(N-p_e-1)\end{array}\right]---\left(1\right)$式中，p_e为所取信号模型的阶数，x(i)表示第i个测量数据，i＝0、1、2、…、N‑2，N为测量数据总个数；C、用残差迭代法计算特征方程系数取信号模型的阶数m为满足p≤m<[N/2]的整数，列出特征方程系数求解方程组，由于特征方程系数求解方程组的方程数大于未知数个数，不能直接求解，形成相应的法方程后采用最小二乘法对未知数进行估计，计算出特征方程的系数a＝[a₁ a₂ … a_m]^T           (2)特征方程系数求解方程组为：$\left[\begin{array}{cccc}x\left(m\right)& x(m-1)& ...& x\left(0\right)\\ x\left(m\mathrm{+1}\right)& x\left(m\right)& ...& x\left(1\right)\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ x(N-1)& x(N-2)& ...& x(N-m-1)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ·\\ ·\\ ·\\ a_m\end{array}\right]=\left[\begin{array}{c}ϵ\left(m\right)\\ ϵ(m+1)\\ ·\\ ·\\ ·\\ ϵ(N-1)\end{array}\right]---\left(3\right)$式(3)对应的法方程为：$\left[\begin{array}{cccc}r(0,0)& r(0,1)& ...& r(0,m)\\ r(1,0)& r(1,1)& ...& r(1,m)\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ r(m,0)& r(m,1)& ...& r(m,m)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ·\\ ·\\ ·\\ a_m\end{array}\right]=\left[\begin{array}{c}{ϵ}_m\\ 0\\ ·\\ ·\\ ·\\ 0\end{array}\right]---\left(4\right)$式中，样本函数r(i,j)的表达式如下：$r(i,j)=\underset{n=m}{\overset{N-1}{\Sigma }}x(n-j)x^{*}(n-i),i=0,1...,m,j=0,1...,m---\left(5\right)$上式中的上标(*)表示共轭；最小误差能量为${ϵ}_p=\underset{j=0}{\overset{m}{\Sigma }}{a}_j\left[\underset{n=m}{\overset{N-1}{\Sigma }}x(n-j)x^{*}\left(n\right)\right]---\left(6\right)$信号模型的阶数为m时的特征方程为：$\underset{i=0}{\overset{m}{\Sigma }}{a}_i{z}^{m-i}=0---\left(7\right)$式中，a₀＝1；为了提高计算精度，采用迭代法求解a；其步骤如下：C1、设迭代计数k＝0，收敛精度δ＝0.00001，最大迭代次数k_m＝30；C2、根据测量数据形成方程组的系数矩阵X和右端向量c；用最小二乘法对特征方程系数求解方程组(3)求解，就是把式(3)写成下式所示的不相容方程组，并令误差平方和最小，$\left[\begin{array}{cccc}x\left(m\right)& x(m-1)& ...& x\left(0\right)\\ x\left(m\mathrm{+1}\right)& x\left(m\right)& ...& x\left(1\right)\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ x(N-1)& x(N-2)& ...& x(N-m-1)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ·\\ ·\\ ·\\ a_m\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ·\\ ·\\ ·\\ 0\end{array}\right]---\left(8\right)$X和c分别为式(8)经过变化后得到的方程组的简写形式的系数矩阵和右端向量：Xa＝c        (9)式中，$X=\left[\begin{array}{cccc}x(m-1)& x(m-2)& ...& x\left(0\right)\\ x\left(m\right)& x(m-1)& ...& x\left(1\right)\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ x(N-2)& x(N-3)& ...& x(N-m-1)\end{array}\right]---\left(10\right)$a＝[a₁ a₂ … a_m]^T           (11)c＝‑[x(m) x(m+1) … x(N‑1)]^T         (12)C3、用高斯消去法求解式下式，得到a的初值a⁽⁰⁾；X^TXa＝X^Tc            (13)C4、把a^(k)代入式(c′)^(k)＝Xa^(k)，得到(c′)^(k)；C5、按照式Δc^(k)＝c‑(c′)^(k)，计算残差Δc^(k)；C6、用高斯消去法求解式X^TXΔa^(k)＝X^TΔc^(k)，得到Δa^(k)；C7、求${\mathrm{\Delta a}}_{max}=\underset{1\le i\le m}{max}\left|{\mathrm{\Delta a}}_i^{\left(k\right)}\right|;$C8、判断Δa_max是否小于δ，如果Δa_max小于δ，转至步骤D；C9、令a^(k+1)＝a^(k)+Δa^(k)；C10、令k＝k+1，如果k<k_m，转至步骤C4；D、计算频率和衰减因子特征方程的系数a求出后，代入特征方程(7)，得到z后，并计算出信号分量的频率和衰减因子；特征方程(7)是一元高次方程，其解为实数或共轭复数，采用QR分解法求解；信号分量的频率和衰减因子计算公式如下：$\left\{\begin{array}{c}{f}_i=arctan[\mathrm{Im}\left(z_i\right)/\mathrm{Re}\left(z_i\right)]/\left(2\pi \Delta t\right)\\ {\alpha }_i=ln\left|z_i\right|/\Delta t\end{array}\right.,i=1,...,m---\left(14\right)$式中，arctan、ln、Im、Re分别为反正切函数、自然对数函数、取复数虚部函数、取复数实部函数；E、计算幅值和初相把求出的z代入以下公式(15)，并根据该式计算出b后，再计算信号分量的幅值和初相；式(15)的方程数也大于未知数个数，不能直接求解，形成相应的法方程，采用最小二乘法对未知数进行估计,计算出信号分量按幅值由大到小排序，取前p个信号分量作为最终的信号分量结果；计算b的公式为：Zb＝x              (15)式中$Z=\left[\begin{array}{cccc}1& 1& ...& 1\\ z_1& z_2& ...& z_m\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ z_1^{N-1}& z_2^{N-1}& ...& z_m^{N-1}\end{array}\right]---\left(16\right)$b＝[b₁ b₂ … b_m]^T            (17)x＝[x(0) x(1) … x(N‑1)]^T         (18)相应的法方程为：Z^TZb＝Z^Tx幅值和初相计算公式为$\left\{\begin{array}{c}{A}_i=\left|b_i\right|\\ {\theta }_i=arctan[\mathrm{Im}\left(b_i\right)/\mathrm{Re}\left(b_i\right)]\end{array}\right.,i=1,...,m---\left(19\right).$