exprony_ma.m

基于prony算法的低频振荡模式识别

%调用prony算法函数
function[F,D,A,theta]=exprony_ma(y,p,sf);
nm=2*p;
n=fix(length(y)/2);
dt=1/sf;%由采样频率计算时间间隔
L=length(y);
M=L/2;
for k=1:nm
x1(:,k)=y(k:M-1+k);
end
for k=1:M
x2(k,:)=-y(nm+k);
end
B=x1\x2;%最小二乘法求解prony多项式系数
B(nm+1)=1;
B1=B(nm+1:-1:1);
V=roots(B1);%求出特征方程的根
F1=abs(log(V))/(2*pi*dt);%初步计算频率
D1=log(abs(V))/dt;%初步计算衰减因子
for k=0:(2*n-1)
Va(k+1,:)=conj(V').^k;%生成范得蒙矩阵Va
end
S1=2*(inv(conj(Va')*Va)*conj(Va')*y');
[F2,I]=sort(F1);%将频率F1从小到大排序
m=0;
for k=1:nm-1
if F2(k)~=F2(k+1)%消除其中的共轭项
continue;
end
m=m+1;
ii=I(k);%提取有效的元素序数
F(m)=F1(ii);
D(m)=D1(ii);
A(m)=abs(S1(ii));
theta(m)=angle(S1(ii));
end