基于prony算法的参数辨识算法的仿真——简化版

基于prony算法的参数辨识算法的仿真的详细版欢迎订阅本博：

https://blog.csdn.net/ccsss22/article/details/115358232

1.问题描述：

 建立如下被测信号：

被测信号中包含四个振荡模态，在数据窗宽度同样为10s的前提下，利用不同的采样频率做普罗尼计算。 

2.部分程序：
 
function X = func_Prony(Signal,dt);

s     = Signal; 
L     = length(s(1:length(Signal))); 
Order = ceil(L/2); 
R     = []; 
K1    = 0;   
K2    = 0; 

%扩展矩阵 
while K1 <= Order 
    K2  = 1; 
    Re = []; 
    while K2 <= Order 
        u  = Order - K2 +1; 
        v  = L - K2; 
        m  = Order - K1 +1; 
        l  = L - K1; 
        r  = sum(s(u:v).*conj(s(m:l))); 
        Re = [Re,r]; 
        K2  = K2 + 1; 
    end 
    R = [R,Re']; 
    K1 = K1 + 1; 
end 

%计算阶数 
Order = func_Order(R); 

%计算相关参数
K2 = Order-1; 
Re = R(2:end,2:K2+1); 
b  = R(:,1); 
b  = b(2:end); 
a  = pinv(Re)*(-b); 

R1 = R(1,:); 
R1 = R1(1:K2+1); 
a1 = [1 a']; 

Ep = sum(R1.*a1); 
P  = [1 a']; 
z  = roots(P); 

%估计序列X 
ks    = 1:K2; 
X(ks) = s(ks); 

for Order = K2+1 : L 
    Lij      = 1:K2; 
    X(Order) = sum(-a'.*s(Order-Lij)); 
    Order    = Order+1; 
end 

Zh = []; 
for mm = 0:L-1 
    Zh = [Zh,z.^mm]; 
end 

Z    = Zh'; 
Z    = conj(Z); 
Zhh  = Z'; 
b    = (Zhh*Z)^(-1)*Zhh*X'; 
%最后得到的四个参数值
A     = abs(b)
f     = angle(z)/2/pi/0.001
a     = log(abs(z))/dt
theta = angle(b)/2/pi/dt

 

3.仿真结论：
       注意，这里论文中你所给的那个公式，貌似有点小错误，这里我们使用了两组公式进行计算，一组是你所提供的公式，一组是我们给的测试数据。

       仿真结果如下所示：

A-27-6