prony1.m

来自「可以用来识别包装系统的总体参数以及原响应函数的初始相位和振幅」· M 代码 · 共 18 行

t=0:0.001:1;
x=5*exp(-0.008*2*pi*50*t).*sin(2*pi*50*t+pi/3.0);
for p=8:2:50
    for i=1:1001-p
        for j=1:p
            q(i,j)=x(p-j+i);
        end
    end
    b=-x(p+1:1001);
    b=b';
    p
    a=pinv(q)*b;
    plot(roots([1,a']),'+');
    hold on
    clear q
    roots([1,a'])
end