关于prony算法

xian2006

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

  X8 Q) r4 V% g9 _所以在这里诚心向各位请教
. a# X: I$ W: A: J- T9 i1 ^; p

1. %% 数据准备
   
2. clear;
   8 l0 V8 I/ ?/ [% h3 T$ x) `
3. clc;
   
4. format long
   / r* j! W9 X1 N( _$ p+ M8 r
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   ; ^* O+ |/ O  p$ q5 P: e7 _
8. f1 = 49;
   ! M  i$ Q+ B" O
9. f2 = 51;
   & N" ]1 y, _! n  c. h
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);
    ! O0 N" _; i8 j% q
14. pe = floor(N/2);
    
15. 
    ' M- O3 v! X/ G: C/ U
16. %% 构造样本矩阵
    . @1 V- S9 ]$ ?" R
17. 
    
18. Re=zeros(pe+1,pe+1);
    # i. \! S* m" b1 Q
19. 
    
20. for i = 2:pe+1
    ; |( r0 p4 F! D- @7 M/ E
21. for j = 1:pe+1
    , b7 A; t( x) _
22. for n = pe:N-1
    5 l# j, p9 r8 r) K+ n. u
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    ' ]( h4 V4 q: r8 M8 }+ }4 s* @  D
24. end
    
25. end
    7 N: P( m. ?9 Y8 }0 U9 G
26. end
    7 R+ V8 a& ]- e, O+ U0 j0 F
27. 
    
28. Re(1,:) = [];
    
29. 
    
30. %% SVD_TLS确定介数p及a
    3 }2 Q4 G  ]" l1 `, a; z
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    9 t9 O) ]6 @, E$ u' t, U0 o- \* D9 T
33. 
    
34. % 求p值
    
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    
37. for i = 1:pe
    + M: M- O- R: y9 h
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    
40. 
    
41. % 归一化比值Ak/A 求p值
    % k6 j* ?1 ?5 }+ Q- Y" h1 Z0 Q1 A
42. sum_k = 0;
    7 U" \, k3 O1 n% X9 B/ k1 ~
43. k = 1;
    
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    % x! K# @  R! @* [. l7 P
45. sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
    , H7 S% _+ D9 s3 s7 @2 E& ?6 m
46. k = k+1;
    
47. end
    
48. p = k-1;
    3 u& C! B; a3 `! J7 |7 T
49. 
    
50. % 求Sp部分
    & m% E7 F2 I* p& M- M7 g
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
    ( Q; \5 O- t7 g, ]* b' d
56. end
    
57. 
    
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    - L0 R. `+ x4 \# d" T! \
60. % SS(i,i) = S(i,i);
    
61. % end
    + B' o( ]/ t' i% Z/ ^
62. % B = U*SS*V;
    
63. 
    % ~4 |. G6 d, F, _
64. % 求Sp逆矩阵
    4 u+ |9 G  G0 W- D& k! I, }
65. inv_Sp=inv(Sp);
    1 L5 Y; q' P4 @9 j$ D  m
66. if isinf(inv_Sp(1,1)) == 1
    % w" K3 [- x* p: p' g" F6 Y
67. inv_Sp = pinv(Sp);
    
68. end
    
69. 
    , B2 s" ?; K: Y$ F
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    
72. 
    
73. %% 求z
    
74. y=[1 a'];
    7 g* O3 H# S* C" a! _* \
75. z=roots(y);
    6 ^. F% R4 P9 B4 m
76. 
    ( M. V: o8 H) k0 Q. c' Y
77. %% 求x的近似值x_j
    
78. %求前p近似值等于测量值 x_j(1:p)
    , g& r  D; f2 i( g* B8 I& l
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
    & g: }1 I8 U( t% w9 H
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    . t9 J' P6 [" E7 W- t
88. end
    + z* e& I2 y( W$ t3 u' K  r
89. end
    8 ~5 {4 V7 [$ Y5 c: I" V: C
90. 
    
91. %% 画图 x、x_j
      d$ B- M- {" L# q2 u+ u6 L6 C$ z
92. hold on;
    7 f5 [  U$ P( ]+ y0 r
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    
95. hold off;
    3 y8 G9 y2 v; P* Y- D4 y$ w/ q1 V
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    : [! F/ C% Z4 [" V
98. 
    * |- s" s* O; C7 b) v8 @* ~: e
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     + b% a9 `3 y" a
101. for i=1:N
     , }. Z8 [+ C+ R' j8 i$ E' ^( B
102. Z(i,:)=z'.^(i-1);
     5 c0 u# i" S; y8 I7 M( |6 @
103. end
     
104. 
     . G3 p+ h0 ^$ T+ Z6 C1 o5 I
105. %求取H
     8 Y9 `% e2 W: V+ u, F
106. H=zeros(p,p);
     8 I! ^% E$ Q2 p9 y: V/ H
107. for i=1:p
     
108. for j=1:p
     & K; n. n* Y5 V* f' R6 N5 f2 N+ E
109. m=(conj(z(i))*z(j));
     
110. H(i,j)=(m^N-1)/(m-1);
     6 V# a" P" p- |( K* K# \% h
111. end
     1 ?0 ~. Y5 x4 V! T$ [$ `
112. end
     
113. 
     1 R. @6 r) y# m: o8 u
114. % 求取b
     
115. b=inv(H)*Z'*x';
     
116. 
     
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     % g3 t4 b% J) @2 }6 m/ `
118. for i = 1:p
     ' Q2 P( Y6 s" f  t* {' ]
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;
     5 J6 y! w; i6 ^6 G& `5 }( l
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323

5 m1 e( M" a1 \2 e, Q2 Y
) g( l6 t% r; X1 z: N, h    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
. q4 t% V" g) z* x: T9 S
. N& F6 u9 z% s2 ^1 _. g
    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
- R" k* P8 n' c% a' R) E& v
$ z- Z- [8 h! o0 G! u
7 t$ E$ i2 k2 K& z    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶