关于prony算法

xian2006

最近一直在搞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 ，增加数据窗长度还是失败，而且有数据稳定性也有问题。
$ F. Z5 a: j+ v1 x$ m2 b但是用mathwork网站上下的pronytool工具箱还是能比较容易的分离出来的，于是我查看pronytool的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
6 K/ w3 Q1 l& q
所以在这里诚心向各位请教
9 u, P: j$ [5 d. W" }( S' ]* E

1. %% 数据准备
   . Y5 X" c; d) S8 @, e) g
2. clear;
   
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   0 j$ u7 ^) H" h# i
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   
8. f1 = 49;
   
9. f2 = 51;
   $ ?) t& p8 I5 R. p1 S, q5 j4 \4 c
10. t=0.0001:0.0001:0.01;
    : ^2 p6 s, t6 ]& W6 E+ k
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;
    / q! M/ n: d: P1 x9 _& Y. U2 Q: s
13. N = length(x);
    & ^2 v! v$ Z( {  M/ E2 q" d4 i0 f
14. pe = floor(N/2);
    
15. 
    
16. %% 构造样本矩阵
    3 u. F0 M6 U* V0 B5 z
17. 
    # t" R4 n4 r& ]& {0 J
18. Re=zeros(pe+1,pe+1);
    
19. 
    ( _4 q9 I2 R. O9 a7 A
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);
    4 [# ^# D6 B( D" G2 N$ C" t4 C
24. end
    5 R, b- b3 ?+ O$ I5 `; ?. ~1 |
25. end
    
26. end
    & [3 _( j; t  t1 q& f
27. 
    
28. Re(1,:) = [];
    & M0 v. T3 e+ }- w5 w: J1 T/ T
29. 
    9 A& u, C! q9 a
30. %% SVD_TLS确定介数p及a
    
31. 
    0 B/ ~. Y+ h" M7 Y
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    3 q) m, B; i9 m2 c
33. 
    2 z0 V2 y' _' C
34. % 求p值
    
35. % 计算全部奇异值平方和
    & @" u/ C) R: l. ]- m
36. sum_all = 0;
    
37. for i = 1:pe
    1 U  a& N; P- F2 ^  v$ R& l
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    # B' `. O5 j( q" |) g% B* D
40. 
    ( y/ R5 }8 B& {& l( \
41. % 归一化比值Ak/A 求p值
    $ m* G/ C1 \, A' X; p
42. sum_k = 0;
    
43. k = 1;
    
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    6 a' h4 Z1 g  E' A" B: ]/ ?/ o
45. sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
    ) M6 r  f- }! {2 i$ a  D
46. k = k+1;
    1 j! @5 n6 o0 A1 u
47. end
    % J# @' ^* _: S0 c
48. p = k-1;
    
49. 
    5 A* T: E0 A1 n) i9 s( U$ I
50. % 求Sp部分
    ( j: q' b) ]6 I# ^
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)';
    3 w7 {' D/ f. ]# u/ D
55. end
    # O( \1 M! I9 ^3 b9 D5 d) t, K
56. end
    
57. 
    7 K% P% s: F! V5 r" L
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    
61. % end
    7 l/ b4 |3 h. R- r* N3 `
62. % B = U*SS*V;
    
63. 
    
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
    ( @' e- Z3 L/ F. R# |, H9 w
66. if isinf(inv_Sp(1,1)) == 1
    ) Q0 _5 E9 @1 ?2 L/ ~) Y$ ^+ b
67. inv_Sp = pinv(Sp);
    
68. end
    
69. 
    
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    6 e' v4 O+ v( W* F$ j$ U9 M
72. 
    
73. %% 求z
    3 s) q. S6 ]+ z) Y. |2 m
74. y=[1 a'];
    9 g! k5 U# [- S! E" j3 Y' Q
75. z=roots(y);
    $ x  z! w( ?' b" \" _6 A: K" ?
76. 
    
77. %% 求x的近似值x_j
    * e7 K0 S0 X% V8 Q% U: d. r
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    - R3 W% t* l. Y; `
80. for i = 1:p
    
81. x_j(i)=x(i);
    
82. end
    # r' R" s" i* O+ _! V
83. 
    
84. %求x的N-p+1个近似值 x_j(p+1:N)
      I5 b0 G# A- d+ A, ], Y2 w3 e
85. for n = p+1:N
    5 _- s4 r/ g; h& E) J% Q
86. for i = 1:p
    . R* t* X! Y1 |# Z# P
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    % ^/ R! f2 R8 l
88. end
    
89. end
    
90. 
    0 A" Q' t2 t0 t* k0 y3 r
91. %% 画图 x、x_j
    - L: @. _- N( U. X; u) y. y  z
92. hold on;
    ! @  a, K# c: Z/ m( K
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
    " r- {3 Q) N8 [1 p/ O
100. Z=zeros(N,p);
     
101. for i=1:N
     
102. Z(i,:)=z'.^(i-1);
     & B. [2 p+ f# Z3 \, B# {
103. end
     
104. 
     " G) l" k; H& m2 q" d6 D
105. %求取H
     
106. H=zeros(p,p);
     
107. for i=1:p
     
108. for j=1:p
     / r# d6 ]2 ^, {0 J
109. m=(conj(z(i))*z(j));
     
110. H(i,j)=(m^N-1)/(m-1);
     : B" \. l4 t- v  ]. U  V2 G
111. end
     
112. end
     # G1 @# `% R4 p- H6 _! B
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     7 t+ r/ M2 P+ m) P  G/ d: |
116. 
     * v7 m7 n8 p" o! C4 [4 v. f" t
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     2 P9 y# M* z0 L! k: `
118. for i = 1:p
     $ V, j4 `" M* r3 h! K
119. Amp(i) = abs(b(i));
     # `5 l  V5 W7 d
120. Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     7 z% z! E0 {- B5 W) s
121. Damp(i) = log(abs(z(i)))*dt;
     ) v% g- f# Y  u7 E' o4 E
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     * S) ?* ^9 t/ {8 B% [7 Q) Y, O. T# n
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
& G/ h. t: p$ _& F- m

% W. |+ H; e5 f6 p    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
0 x+ ]: s, P* k( M
+ O6 j  H* u8 ~4 I5 s6 D1 d
( r. B8 x  z  A& n. \    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
/ j4 t" g: [5 q+ `) w& Z4 `

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

nanshugong

我帮学长顶顶