limprony.m

来自「功率谱估计程序」· M 代码 · 共 30 行

%prony算法，选用电力系统谐波信号作为输入信号。参考书目：现代信号处理(第二版），张贤达，清华大学出版社。122－124页
clc
clear
x=0:199;
%低频振荡信号
x=220*exp(-0.25*x*0.05).*cos(2*pi*1.5*x*0.05+pi)+110*exp(-0.35*x*0.05).*cos(2*pi*1.0*x*0.05+pi/2);p=4;N=200;ts=0.05;
for i=0:p
    for j=0:p
        s=0;
        for n=p:N-1
            s=s+x(n-j+1)*x(n-i+1);
        end
        r(i+1,j+1)=s;
    end
end
r1=r(2:end,2:end);
r2=r(2:end,1);
a=-inv(r1)*r2;
a=[1,a'];
z=roots(a);
z=z';
for n=0:N-1
    z1(n+1,:)=z.^n;
end
z=z';
b=(z1'*z1)^(-1)*z1'*x';
    A=abs(b)
    s=angle(b)*180/pi
    erfa=log(abs(z))/ts
  f  = angle(z)/2/pi/ts