关于prony算法

最近一直在搞prony算法，资料是张贤达的《现代信号处理》和束红春的《电力工程信号处理应用》，这两本书都有关于prony扩展算法的内容，甚至有相关步骤，我也根据这些内容自己编了程序。为了测试 在信号x=160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
取时窗10ms 100个数据点，但是分离不出这两信号，后面尽管修改程序自己设定介数p ，增加数据窗长度还是失败，而且有数据稳定性也有问题。
但是用mathwork网站上下的pronytool工具箱还是能比较容易的分离出来的，于是我查看pronytool的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
所以诚心向各位请教。

1. %% 数据准备
   
2. clear;
   
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   
7. % t = t(201:500); exp_matrix.m (2.04 KB, 下载次数: 22935)
   2011-11-15 22:15 上传

   点击文件名下载附件
   程序
   
8. f1 = 49;
   
9. f2 = 51;
   
10. t=0.0001:0.0001:0.01;
    
11. x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
    
12. dt = 0.0001;
    
13. N = length(x);
    
14. pe = floor(N/2);
    
15. 
    
16. %% 构造样本矩阵
    
17. 
    
18. Re=zeros(pe+1,pe+1);
    
19. 
    
20. for i = 2:pe+1
    
21.     for j = 1:pe+1
    
22.         for n = pe:N-1
    
23.         Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    
24.         end
    
25.     end
    
26. end
    
27. 
    
28. Re(1,:) = [];
    
29. 
    
30. %% SVD_TLS确定介数p及a
    
31. 
    
32. [U,S,V] = svd(Re);       %%%%%%奇异值分解
    
33. 
    
34. % 求p值
    
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    
37. for i = 1:pe
    
38.     sum_all = sum_all+S(i,i)^2;
    
39. end
    
40. 
    
41. % 归一化比值Ak/A 求p值
    
42. sum_k = 0;
    
43. k = 1;
    
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    
45.     sum_k = sum_k+S(k,k)^2;                      %%%%%%%计算k个奇异值平方和
    
46.     k = k+1;
    
47. end
    
48. p = k-1;
    
49. 
    
50. % 求Sp部分
    
51. Sp=zeros(p+1,p+1);                        %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    
52. for j = 1:p
    
53.     for i = 1:(pe+1-p)
    
54.        Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    
55.     end
    
56. end
    
57. 
    
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    
60. %     SS(i,i) = S(i,i);
    
61. % end
    
62. % B = U*SS*V;
    
63. 
    
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
    
66. if isinf(inv_Sp(1,1)) == 1
    
67.     inv_Sp = pinv(Sp);
    
68. end
    
69. 
    
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    
72. 
    
73. %% 求z
    
74. y=[1 a'];
    
75. z=roots(y);
    
76. 
    
77. %% 求x的近似值x_j
    
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    
80. for i = 1:p
    
81.     x_j(i)=x(i);
    
82. end
    
83. 
    
84. %求x的N-p+1个近似值  x_j(p+1:N)
    
85. for n = p+1:N
    
86.     for i = 1:p
    
87.         x_j(n)=x_j(n)-a(i)*x_j(n-i);
    
88.     end
    
89. end
    
90. 
    
91. %% 画图 x、x_j
    
92. hold on;
    
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    
95. hold off;
    
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    
98. 
    
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     
101. for i=1:N
     
102.     Z(i,:)=z'.^(i-1);
     
103. end
     
104. 
     
105. %求取H
     
106. H=zeros(p,p);
     
107. for i=1:p
     
108.     for j=1:p
     
109.         m=(conj(z(i))*z(j));
     
110.         H(i,j)=(m^N-1)/(m-1);
     
111.     end
     
112. end
     
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     
116. 
     
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     
119.     Amp(i) = abs(b(i));
     
120.     Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     
121.     Damp(i) = log(abs(z(i)))*dt;
     
122.     Pha(i) = atan(imag(b(i)/real(b(i))));
     
123. end

复制代码

是不是算法程序出了问题 求帮忙

把错误贴出来以便会员解答

可以运行，就是打不到要的效果，分离不出两个频率相近的信号

可以共享一下mathwork网站上下的pronytool工具箱么644682383@qq.com，真心感谢。一同学习。。

大家有prony工具箱的中文说明吗，我分析的时候频率是成对出现的，应该怎样处理啊

yyl_liang 发表于 2013-3-28 10:35
大家有prony工具箱的中文说明吗，我分析的时候频率是成对出现的，应该怎样处理啊
...

遇到了相同的问题，谁有工具箱的中文说明啊

我在mathworks网站上下了pronytool工具箱，结果用半天不会用 求大神指导下我啊  看英文看不大懂 大神加我qq：765481889

同求工具箱的中文说明

yyl_liang 发表于 2013-3-28 10:35
大家有prony工具箱的中文说明吗，我分析的时候频率是成对出现的，应该怎样处理啊
...

请问你的问题是如何解决的啊？

yyl_liang 发表于 2013-3-28 10:35
大家有prony工具箱的中文说明吗，我分析的时候频率是成对出现的，应该怎样处理啊
...

频率成对出现的问题，而且相位等也是成对的，正负号不同，怎么解释啊

楼主你好！我检验了你的Re发现结果第一行全是0，请问为什么？谢谢

prony工具箱怎么使用呢？有中文说明吗？

同样再弄prony,遇到了同样问题

精度上还是差很多啊。怎么提高精度啊。你这个算法流程和步骤应该是没有问题的。

huanghun0912 发表于 2013-1-14 10:28
可以共享一下mathwork网站上下的pronytool工具箱么，真心感谢。一同学习。。

同求  487404058@qq.com  万分感谢

楼主这个程序是工具箱上的，还是书上的呢？

楼主可不可以分享下下载的Prony工具箱安装包呢？万分感谢！我邮箱：2535713512@qq.com

丽红 发表于 2018-11-29 11:40
同求    万分感谢

你好，请问你收到了PRONY的安装包了吗？可以分享下咩？感谢~我邮箱：2535713512@qq.com

你好，请问一下点击Start运行时有没有出现如下bug，如果有的话怎样解决，谢谢！
pronytool
Undefined function or variable 'xlate'.

Error in pronytoolbar (line 24)
mlabel = get(h1(j),xlate('TooltipString'));

Error in pronyanalysistool>pronyanalysistool_OpeningFcn (line 62)
pronytoolbar(hObject,'on');

Error in gui_mainfcn (line 220)
feval(gui_State.gui_OpeningFcn, gui_hFigure, [], guidata(gui_hFigure), varargin{:});

Error in pronyanalysistool (line 42)，

john.bham 发表于 2019-6-3 02:01
你好，请问一下点击Start运行时有没有出现如下bug，如果有的话怎样解决，谢谢！
pronytool
Undefined func ...

我也出现这个错误了，请问你解决了吗

大佬，这个PRONY的输入参数是采样时间和函数吗？

x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
dt = 0.0001;
N = length(x);
pe = floor(N/2);
这个PRONY的输入分别是x、dt、N、pe这四个参数吗？我有仿真数据，但是不知道怎么拿PRONY来分析我我的仿真数据，x代表要分析信号的函数、dt代表采样时间、N是代表数据窗的长度、pe是代表取整。如何输入这四个参数对我的仿真数据进行分析呢？我的仿真数据类型是exccel中的数据。

%% 数据准备
clear;
clc;
format long
% load('1000kV示范工程线路','t','vX0043a')
% x = vX0043a(201:500);
% t = t(201:500);

f1 = 40;
f2 = 60;
t=0.001:0.001:0.05;
x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
dt = 0.001;
N = length(x);
pe = floor(N/2);

%% 构造样本矩阵

Re=zeros(pe+1,pe+1);

for i = 2:pe+1
    for j = 1:pe+1
        for n = pe:N-1
        Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
        end
    end
end

Re(1,:) = [];

%% SVD_TLS确定介数p及a

[U,S,V] = svd(Re);       %%%%%%奇异值分解

% 求p值
% 计算全部奇异值平方和
sum_all = 0;
for i = 1:pe
    sum_all = sum_all+S(i,i)^2;
end

% 归一化比值Ak/A 求p值
sum_k = 0;
k = 1;
while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    sum_k = sum_k+S(k,k)^2;                      %%%%%%%计算k个奇异值平方和
    k = k+1;
end
p = k-1;

% 求Sp部分
Sp=zeros(p+1,p+1);                        %%%%%%s生成（p+1）X（p+1）维矩阵Sp
for j = 1:p
    for i = 1:(pe+1-p)
       Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    end
end

% SS = zeros(pe,pe+1);
% for i = 1:p
%     SS(i,i) = S(i,i);
% end
% B = U*SS*V;

% 求Sp逆矩阵
inv_Sp=inv(Sp);
if isinf(inv_Sp(1,1)) == 1
    inv_Sp = pinv(Sp);
end

% 求a
a=inv_Sp(2:p+1,1)/inv_Sp(1,1);

%% 求z
y=[1 a'];
z=roots(y);

%% 求x的近似值x_j
%求前p近似值等于测量值 x_j(1:p)
x_j=zeros(N,1);
for i = 1:p
    x_j(i)=x(i);
end

%求x的N-p+1个近似值  x_j(p+1:N)
for n = p+1:N
    for i = 1:p
        x_j(n)=x_j(n)-a(i)*x_j(n-i);
    end
end

%% 画图 x、x_j
hold on;
plot(t,x,'k');
plot(t,x_j,'r');
hold off;

%% 求取 b=inv(H)*Z'*x_j

% 求取N X p维vandermode矩阵Z
Z=zeros(N,p);
for i=1:N
    Z(i,:)=z'.^(i-1);
end

%求取H
H=zeros(p,p);
for i=1:p
    for j=1:p
        m=(conj(z(i))*z(j));
        H(i,j)=(m^N-1)/(m-1);
    end
end

% 求取b
b=inv(H)*Z'*x';

%% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
for i = 1:p
    Amp(i) = abs(b(i))
    Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt)
    Damp(i) = log(abs(z(i)))*dt
    Pha(i) = atan(imag(b(i)/real(b(i))))
end



Amp =

  75.435264220318473  75.435264220318459  79.821614453380022  79.821614453380022


Fre =

  60.000000000644370 -60.000000000644370  39.999999999685848 -39.999999999685848


Damp =

   1.0e-05 *

  -0.299999999623960  -0.299999999623960   0.000000000270939   0.000000000270939


Pha =

   0.475722414206541  -0.475722414206541   0.654158664114685  -0.654158664114685

鄙人觉得您时间应该从零开始