function?[jxSNRresult]=myprony（xpdt）
%?實現《電網諧波與間諧波檢測的分段Prony算法》郭成?1
%?實現《現代信號處理》張賢達P122-124
%?實現《Prony算法在諧波、間諧波參數辨識中的應用》式（12）?楊玉坤
%?intput?x:原始數據
%?intput?p:prony模型階數
%?intput?dt:采樣間隔
%?output?jx:擬合數據
%?output?Am：振幅
%?output?Qm：相位
%?output?am：衰減因子
%?output?fm：振蕩頻率
%?output?result=[AmQmamfm]

N=length（x）;

%%?廣義求逆解ai
%建立矩陣X
X=zeros（N-pp+1）;
for?i=0:N-p-1
????X（i+1:）=x（p+i+1:-1:i+1）;
end
[mn]=size（X）;
Xa=X（:2:n）;
Xb=-X（:1）;
%最小二乘法求ai
ai=pinv（Xa）*Xb;
ai=ai‘;

%%?解特征方程，得zi
zi=roots（[1ai]）;
zi=zi‘;

%%?最小二乘法求bm
%建立矩陣Z
Z=ones（Np）;
for?i=0:N-1
????Z（i+1:）=zi.^i;
end
%最小二乘法求bm
bm=pinv（Z）*x‘;

%%?求其他參數
%賦值Am
Am=abs（bm）;
%頻率fm及幅角修正
fm=atan（imag（zi）./real（zi））;
fm=fm‘;
for?i=1:p
????if?real