关于prony算法

xian2006

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

5 Q3 t/ I3 O  f. e所以在这里诚心向各位请教
8 y3 B- d& j7 ?: Y  t" [" O- S

1. %% 数据准备
   0 C! s" ^, c; j2 K9 k$ j+ e* Y9 s7 N
2. clear;
   
3. clc;
   
4. format long
     Y; E& o! S0 V, q  u4 Y" ?
5. % load('1000kV示范工程线路','t','vX0043a')
   7 y( [- p6 r! w
6. % x = vX0043a(201:500);
   / F: o; z- w1 E1 `
7. % t = t(201:500);
   + i- B" |0 j7 Q1 c2 M
8. f1 = 49;
   9 s/ ?8 R6 x/ s* x6 ?4 y
9. f2 = 51;
   ' Z! A* z' m2 W1 M" ]! }+ U, E
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);
    8 {1 k+ z# Q- R
14. pe = floor(N/2);
    & U. ^# f) i. [! m, `
15. 
    
16. %% 构造样本矩阵
    # {1 A) ?' O; D2 \
17. 
    
18. Re=zeros(pe+1,pe+1);
    ; h1 t! d& @2 G  ^% i3 k2 D2 r
19. 
    0 H/ \  o6 i/ W
20. for i = 2:pe+1
    
21. for j = 1:pe+1
    ; ]9 q; L+ V9 S
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
    " R. c% i( ^3 r" \' i) @
27. 
    - y3 {" Y+ d* v$ Q1 h7 `
28. Re(1,:) = [];
    ' a  C; a5 _" x1 K
29. 
    3 i3 H( y7 b' M" A; X
30. %% SVD_TLS确定介数p及a
    
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    
34. % 求p值
    5 y7 c  o& r, Y- C5 N
35. % 计算全部奇异值平方和
    $ L6 l& U. K1 E  Z1 e; C
36. sum_all = 0;
    & K3 i/ e' J" V" W/ v
37. for i = 1:pe
    " G! \% r! _$ K5 O% @5 Y( U
38. sum_all = sum_all+S(i,i)^2;
    5 {$ p( a4 i6 Z' P6 I) t" Z
39. end
    ! Q* r% g+ C; u6 R; d
40. 
    
41. % 归一化比值Ak/A 求p值
    ( g; h2 F$ D8 S  O3 G. i$ m8 {3 ]
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
    . f8 [. a/ q+ r. n
48. p = k-1;
    
49. 
    0 H* h# k" |: M/ P/ ^. X
50. % 求Sp部分
    . }) y- _3 U$ b/ e2 O7 T" G* c
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    9 j- h+ v  y( \. a+ d
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. 
    $ h: w7 j: ]; ~7 D9 T+ c1 V; |
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逆矩阵
    . B+ U6 D4 {2 ?) j& g
65. inv_Sp=inv(Sp);
    
66. if isinf(inv_Sp(1,1)) == 1
    
67. inv_Sp = pinv(Sp);
    # B* _$ X$ D* A  E
68. end
    ( K& E" O: |& h: b) _, k) T# X& X
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. 
    , W( u9 V( N9 M! d  P2 M, S
77. %% 求x的近似值x_j
    & j3 ^. H- r* ]; Q/ d! I
78. %求前p近似值等于测量值 x_j(1:p)
    ' R) B8 z' D$ h( }: J6 [& `
79. x_j=zeros(N,1);
    0 o, H6 z4 b+ }6 g7 }$ u8 r
80. for i = 1:p
    8 o( S& r& T2 G, w# u6 h
81. x_j(i)=x(i);
    
82. end
      O- W  w. Z/ B+ a; M7 v
83. 
    ' V& }, J5 W4 i2 t" x
84. %求x的N-p+1个近似值 x_j(p+1:N)
    6 ]( r8 Y% y6 D6 L" ?; ~
85. for n = p+1:N
    2 t* [3 K4 T& ?. N& @9 z
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    3 M! g; }3 M" T8 U' p2 T6 a
88. end
    
89. end
    
90. 
    : z, x; e' j9 W6 R9 m! i8 k
91. %% 画图 x、x_j
    
92. hold on;
    
93. plot(t,x,'k');
    9 K* P8 J5 A# [
94. plot(t,x_j,'r');
    
95. hold off;
    . J9 \& K5 l$ [& i5 O$ h0 I
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    
98. 
    $ P! d1 _, U% E, W& k  L6 ^
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     
101. for i=1:N
     
102. Z(i,:)=z'.^(i-1);
     
103. end
     & f/ G, F3 E% G) K
104. 
     . U5 B4 N* `' w2 X2 |+ B6 W
105. %求取H
     " O5 [6 W$ a- N4 r4 |$ T1 d. \
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
     ) J9 D, m( Q4 V
113. 
     
114. % 求取b
     ! Z9 I7 B' Z, A
115. b=inv(H)*Z'*x';
       K- o# ]$ {2 x9 ?1 H" _
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);
     2 y. d* a  I, l2 ^% d: m# j
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
# B0 F$ g; G, A

1 e1 {8 Z) S) e6 A    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
# W! E! i1 B# }& K
" J" e4 r, A/ A5 v* @4 m# i  h
* b; r# Y2 G: L( C6 }9 D6 V: Z8 Q0 A    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
0 A% R5 e9 y/ Y. `- u
( w* y) U* G: }2 S: ^1 N7 H& Y( y
    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶