prony算法的matlab程序实现

function [FDAtheta]=exprony_ma（ypsf）
x=0:0.05:10;
y=220*exp（-0.25*x）.*cos（2*pi*1.5*x+pi）+110*exp（0.25*x）.*cos（2*pi*1.0*x+pi/2）;
p=6;
sf=1/0.05;
nm=2*p; %%% p为模型的阶数
n=fix（length（y）/2）;
dt=1/sf;                        %%% 由采样频率计算时间间隔
h=y（11:2*n）‘;
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;           %% 最小二乘法求解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）