关于prony算法

xian2006

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

所以在这里诚心向各位请教

1. %% 数据准备
   
2. clear;
   7 z4 }% J" p8 A; ^& _
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   2 A: R# L5 ~1 u6 |8 V/ l; X
8. f1 = 49;
   
9. f2 = 51;
   
10. t=0.0001:0.0001:0.01;
    6 t$ b, M, Q0 j; o
11. x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
    ) n4 m8 v0 a& C2 b! l3 ~3 `$ `* R
12. dt = 0.0001;
    * Z' u, u  ~: L, t
13. N = length(x);
    
14. pe = floor(N/2);
      l9 q, \0 s5 \/ P  Q$ h  _! k* F
15. 
    
16. %% 构造样本矩阵
    0 e! r; h2 N( _0 P" E) e
17. 
    7 M# R( q4 B4 m1 Q
18. Re=zeros(pe+1,pe+1);
    
19. 
    " `* V( S8 ]3 T% H; W4 P
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);
    ! k% u% q; D: H
24. end
    % j6 ^. d3 T& `3 `$ N* d
25. end
    
26. end
    # t- w1 f+ O9 C% b7 h8 l
27. 
    
28. Re(1,:) = [];
    
29. 
    ' ]5 m' O) ]. {* x/ `) O; [
30. %% SVD_TLS确定介数p及a
    
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    
34. % 求p值
    0 i; o4 d, d/ f4 g% c- m# G
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    ' [) l4 w' l7 S1 a2 j
37. for i = 1:pe
      S' T7 v# k+ G% R; \+ l9 ^$ v* e/ r' ?
38. sum_all = sum_all+S(i,i)^2;
    , K7 C7 G# T3 z2 P1 e% n
39. end
    ! p8 T0 @& ^$ i' A
40. 
    3 ?% x4 |2 \* w( @. B  o8 _, B! e9 `
41. % 归一化比值Ak/A 求p值
    * c3 e' X& S1 [4 k3 y
42. sum_k = 0;
    3 k6 m/ z6 s# a6 C
43. k = 1;
    . ^+ n  |  q: t  c9 D" Y
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
    1 H( A* [5 a" r1 B
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
    % r4 Q, a" Z1 P! j5 n8 I
56. end
    
57. 
    4 x$ v+ S4 m7 J
58. % SS = zeros(pe,pe+1);
    $ t, f/ ]' k; W) n
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    / p  E: h: O% V; z2 P7 Y
61. % end
    
62. % B = U*SS*V;
    7 U; u/ j$ t; E# Q9 f& C
63. 
    
64. % 求Sp逆矩阵
    8 J4 N( K# B& ]9 W
65. inv_Sp=inv(Sp);
    4 }5 v9 F. y5 m) H: X( |
66. if isinf(inv_Sp(1,1)) == 1
    ! Y1 @5 {/ D$ a# I" u
67. inv_Sp = pinv(Sp);
    - I# ]7 F+ n: v
68. end
    
69. 
    - c7 o$ K2 |; @" s6 W
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    6 P: g8 d( _( @% V6 |* Z
72. 
    
73. %% 求z
    
74. y=[1 a'];
    
75. z=roots(y);
    9 O0 m0 Z0 w; m" e5 `+ \
76. 
    
77. %% 求x的近似值x_j
    2 Z& x3 v9 S6 A( c: d7 Q" o
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    , `& n* Z3 J1 M
80. for i = 1:p
    
81. x_j(i)=x(i);
    
82. end
    
83. 
    
84. %求x的N-p+1个近似值 x_j(p+1:N)
    8 i! e4 m' ^6 a3 o
85. for n = p+1:N
    3 S0 |$ D! e( {8 O/ u; l
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    7 y) I& }3 {3 [2 V
88. end
    
89. end
    
90. 
    
91. %% 画图 x、x_j
    / ~( i) {, G" F5 q
92. hold on;
    * }+ y3 ]& m) q' j/ K& y' n! N
93. plot(t,x,'k');
    # U) z+ @! t, S/ `
94. plot(t,x_j,'r');
    
95. hold off;
    / B; Z1 e/ |' ~3 `
96. 
    ) u+ a( ~3 I' Y4 g
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. 
     3 @+ N) E; W2 v% t
105. %求取H
     
106. H=zeros(p,p);
     ) c/ c' Y# f; W1 D4 A
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);
     * I1 K& i9 b$ y  V
111. end
     
112. end
     ! x) Y$ _$ @9 V; |: Q5 M9 r, [9 h( E0 w
113. 
     . R5 C  |: e' \0 A- B" w5 B
114. % 求取b
     
115. b=inv(H)*Z'*x';
     6 R' H! B2 n; R
116. 
     + l# o  q9 ]2 W& R
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     + R, ?$ @! R" V
118. for i = 1:p
     ) V1 x4 q, K% K2 s
119. Amp(i) = abs(b(i));
     % n3 _1 t( S' F; o
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

复制代码

李真

楼主加油，我们都看好你哦。

Deallee

这个。。。太复杂。额额。不懂哟

xian2006

有人做过这个的话帮个忙啊

tianlangistl

好专业，看不懂帮顶

sjx0323

我的学长似乎在做这个啊！他可能会，帮你请教一下

sjx0323

我的学长似乎在做这个啊！他可能会，帮你请教一下

xian2006

回复 6# sjx0323


" p# J3 M" ~) X! ?$ {4 c, J  a1 [    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006

9 O2 S( V5 y- d# }
    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
' @( v1 `. g6 E1 Z) `

    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶