Prony算法谐波分析matlab函数

Prony算法谐波分析matlab函数,已测试无误，程序使用方法及定义： function [jx,SNR,result]=myprony(x,p,dt) % 实现《电网谐波与间谐波检测的分段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=[Am,Qm,am,fm] ……

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