exprony_ma.m

Prony方法

function [F,D,A,theta]=exprony_ma(y,p,sf)
% F is frequency, D is decaying factor, A is amplitute, theta is phase
% angle.

nm=2*p;
n=fix(length(y)/2);  %取整
dt=1/sf;
h=y;  %
L=length(h);
M=L/2;
for k=1:nm
    x1(:,k)=h(k:M-1+k);
end
for k=1:M
    x2(k,:)=-h(nm+k);
end

B=x1\x2;
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;
end
% % temp=2*(inv(conj(Va')*Va)*conj(Va'));
% % size(temp)
S1=2*(inv(conj(Va')*Va)*conj(Va')*h);

[F2,I]=sort(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));
    % % % A=0;
    % % % theta=0;
end