Prony算法

普罗尼算法（prony algorithm）简介：

　　传统的信号谱分析法如傅里叶算法等，基本上只能分析稳定信号，而对于动态信号却无能为而立。 为解决这个问题，1795年，Prony提出了用指数函数的一个线性组^[1]合来描述等间距采样数据的数学模型，经过适当扩充，形成了能够估算给定信号 的频率、衰减、幅值的prony算法。

优点：

　　在电力系统信号分析尤其是低频振荡分析中显示出良好的应用前景。Prony算法能从均匀采样信号中提取有价值的信息,可分解出一系列按指数衰减的正弦曲线。该算法需要计算最小二乘法拟合和高次代数方程求复根。当信号中嵌入噪声时,Prony求解过程会因方程病态而失败。如输入数据太长,或预估有效模态项数较多,就会造成计算数字溢出。本文方法是:对太长的数据系列采用分段平均值压缩滤波,在求解"法方程"时,能自动确定有效模态的项数,无需SVD奇异值计算,从而大大简化Prony算法、扩大适应范围并确保解算的质量。Prony算法是用一组指数项的线性组合来拟合x(t)=x(n×△τ)等间距采样数据的方法,可以从中分析出信号的幅值B、相位β、阻尼因子σ、频率f等信息。即x(t)=i=ki=1ΣBieσitcos(2πfit+βi)(1)Prony算法实质是对采样数据求拟合公式,用于曲线插值或外推评估.

算法：

　　把时域信号分解成衰减信号：

 

步骤：a、构造离散线性预测模型 b、寻找线性预测模型的特征根 c、利用步骤b决定每一种模型对应的幅值和初相角。

代码：

%% 数据准备   clear;   clc;   format long   % load('1000kV示范工程线路','t','vX0043a')   % x = vX0043a(201:500);   % t = t(201:500);   f1 = 49;   f2 = 51;   t=0.0001:0.0001:0.01;   x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);   dt = 0.0001;   N = length(x);   pe = floor(N/2);  

%% 构造样本矩阵  

Re=zeros(pe+1,pe+1);  

for i = 2:pe+1       for j = 1:pe+1           for n = pe:N-1           Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);           end       end   end  

Re(1,:) = [];  

%% SVD_TLS确定介数p及a  

[U,S,V] = svd(Re);       %%%%%%奇异值分解  

% 求p值   % 计算全部奇异值平方和   sum_all = 0;   for i = 1:pe       sum_all = sum_all+S(i,i)^2;   end  

% 归一化比值Ak/A 求p值   sum_k = 0;   k = 1;   while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe       sum_k = sum_k+S(k,k)^2;                      %%%%%%%计算k个奇异值平方和       k = k+1;   end   p = k-1;  

% 求Sp部分   Sp=zeros(p+1,p+1);                        %%%%%%s生成（p+1）X（p+1）维矩阵Sp   for j = 1:p       for i = 1:(pe+1-p)          Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';       end   end  

% SS = zeros(pe,pe+1);   % for i = 1:p   %     SS(i,i) = S(i,i);   % end   % B = U*SS*V;  

% 求Sp逆矩阵   inv_Sp=inv(Sp);   if isinf(inv_Sp(1,1)) == 1       inv_Sp = pinv(Sp);   end  

% 求a   a=inv_Sp(2:p+1,1)/inv_Sp(1,1);  

%% 求z   y=[1 a'];   z=roots(y);  

%% 求x的近似值x_j   %求前p近似值等于测量值 x_j(1:p)   x_j=zeros(N,1);   for i = 1:p       x_j(i)=x(i);   end  

%求x的N-p+1个近似值  x_j(p+1:N)   for n = p+1:N       for i = 1:p           x_j(n)=x_j(n)-a(i)*x_j(n-i);       end   end  

%% 画图 x、x_j   hold on;   plot(t,x,'k');   plot(t,x_j,'r');   hold off;  

%% 求取 b=inv(H)*Z'*x_j  

% 求取N X p维vandermode矩阵Z   Z=zeros(N,p);   for i=1:N       Z(i,:)=z'.^(i-1);   end  

%求取H   H=zeros(p,p);   for i=1:p       for j=1:p           m=(conj(z(i))*z(j));           H(i,j)=(m^N-1)/(m-1);       end   end  

% 求取b   b=inv(H)*Z'*x';  

%% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha   for i = 1:p       Amp(i) = abs(b(i));       Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);       Damp(i) = log(abs(z(i)))*dt;       Pha(i) = atan(imag(b(i)/real(b(i))));   end