关于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 ，增加数据窗长度还是失败，而且有数据稳定性也有问题。
但是用mathwork网站上下的pronytool工具箱还是能比较容易的分离出来的，于是我查看pronytool的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
; [0 }' e0 o! A! p; v
( R& s" c- r( _* ^9 ^0 {0 E& w所以在这里诚心向各位请教

1. %% 数据准备
   . N7 v. @( d) [% g6 G
2. clear;
   
3. clc;
   
4. format long
   7 e" I( _0 {4 r
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   * v3 E+ |5 I/ A% X  f8 N
7. % t = t(201:500);
   
8. f1 = 49;
   9 L7 a: J. V: d
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;
    9 D6 ~1 w5 k1 t  ?  P+ A' v
13. N = length(x);
    . X0 ]; }" b9 J+ n/ \" t1 f
14. pe = floor(N/2);
    : q) N( M/ b* n% }) C8 ~0 i  M) g/ X1 E
15. 
    
16. %% 构造样本矩阵
    
17. 
    
18. Re=zeros(pe+1,pe+1);
    
19. 
    - m- B4 E$ z- C8 w- d: y0 a8 N
20. for i = 2:pe+1
    
21. for j = 1:pe+1
    
22. for n = pe:N-1
    $ P) `# h4 c6 l0 H
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    
24. end
    
25. end
    ; Q9 ^$ f1 A. @2 H
26. end
    
27. 
    4 M/ ^/ a# O- c+ }  g' H
28. Re(1,:) = [];
    $ J7 A8 [; B1 f: ~4 r  J$ M' @  y
29. 
    . g# {& X' P+ A9 H
30. %% SVD_TLS确定介数p及a
    
31. 
    - R. B* t" l  n: c% M
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    ) `1 \, A  j. _- X# j8 v
34. % 求p值
    
35. % 计算全部奇异值平方和
    ) x- a4 Z3 ?: y2 g( [' k
36. sum_all = 0;
    0 a, m- w2 X! S- E' o
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    + Z- {' i" }! p- p& u; O4 A! h
39. end
    % ~8 x, i% U- j+ z! r
40. 
    
41. % 归一化比值Ak/A 求p值
    
42. sum_k = 0;
    6 w  \, ?, m: O- [1 ?
43. k = 1;
    ; c3 R' l$ X! q, e9 f; D  w
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;
    ( G4 P3 x8 O. C" m7 o! X
47. end
    ; z$ O6 z$ D# \9 ~
48. p = k-1;
    / A2 t, Z) L& i8 x4 }) R% h2 |
49. 
    
50. % 求Sp部分
    
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    ! w0 f3 q$ k* \1 p9 E& G
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
    4 b6 V& e5 t# l5 Y8 f; @; Q; K
56. end
    * k( o5 ~2 q0 A8 S* L6 V
57. 
    . q# }6 O7 p. @- ^! D$ w: V
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    4 `# Q1 g( ?' `3 b) M  O; F
60. % SS(i,i) = S(i,i);
    . F% w& ?1 B" c  w$ l3 l3 E
61. % end
    / `: w; i. \4 Y' c. t
62. % B = U*SS*V;
    ( t% @; {# e5 }5 H. |2 V
63. 
    $ E; [& T& W" D8 }
64. % 求Sp逆矩阵
    3 I$ g# ]5 K# I9 }6 y/ z
65. inv_Sp=inv(Sp);
    * f+ r. ^: U) O& u
66. if isinf(inv_Sp(1,1)) == 1
    * f5 e# d# b, n0 N
67. inv_Sp = pinv(Sp);
    
68. end
    
69. 
    6 ], c" |. V* R" g" C8 a% v
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    4 X, @3 l( [- l$ X" U
72. 
    8 J2 f( f9 C4 E. o2 j
73. %% 求z
    # a, P3 Q0 O% @' f6 U7 I% r* n- Z
74. y=[1 a'];
    # E; C5 U  \/ V- L
75. z=roots(y);
    
76. 
    2 f7 f2 u/ y$ A, s" `3 S
77. %% 求x的近似值x_j
    * A  ^4 B+ \1 W& h; e; `7 I% i, `4 }
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    ) a! I  }9 k( R; k& y! \" M9 [
80. for i = 1:p
    
81. x_j(i)=x(i);
    
82. end
    
83. 
    0 w+ X9 C+ M, o
84. %求x的N-p+1个近似值 x_j(p+1:N)
    
85. for n = p+1:N
    9 h; F4 ^. s; ?8 j
86. for i = 1:p
    8 y$ ^# y7 d$ U* w0 W; g% ^- c; D
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    
88. end
    1 k+ L, B" v& f5 W, D  X* C7 ]
89. end
    - [0 ^, x9 }, p: D9 b; ^
90. 
    - e6 |' i# q5 K; a  e( u
91. %% 画图 x、x_j
    
92. hold on;
    6 U9 H6 _+ v9 G5 B& T1 C* d2 u
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    . A9 ?  A1 N4 ^( r( n. I: N  j3 o4 s
95. hold off;
    * r$ c' V$ _: `; n+ ^3 m8 S& q2 m
96. 
      N! o6 K2 D5 ]
97. %% 求取 b=inv(H)*Z'*x_j
    ! F$ B  {' H& Y' S; U& H
98. 
    7 u8 ~- M, \- v0 r
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     / a# ~9 W2 o; L5 ?" d
101. for i=1:N
     
102. Z(i,:)=z'.^(i-1);
     
103. end
     , M+ d5 k/ c" r: q0 @1 Q
104. 
     
105. %求取H
     
106. H=zeros(p,p);
     5 U/ ]0 V5 ?; C5 _+ d4 U0 t& W
107. for i=1:p
     
108. for j=1:p
     2 ^( |2 q) I/ Q/ g( r) O- [$ _
109. m=(conj(z(i))*z(j));
     
110. H(i,j)=(m^N-1)/(m-1);
     9 S- C7 V# W9 D, M
111. end
     
112. end
     - c2 g* A. t$ K2 @
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     
116. 
     
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     4 v, v7 }: E  F7 X$ O' W
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;
     , S# L1 f% e" @2 s
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323

1 G5 F7 i0 O7 N4 R
    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006


    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
  e9 ^% {+ l" u6 i

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

nanshugong

我帮学长顶顶