关于prony算法

xian2006

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

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

1. %% 数据准备
   ( S$ U  [  X1 [- L% Q% Y6 L- B
2. clear;
   8 k0 A& o0 h/ @' o$ v1 x! N
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   $ T. C' q; F: H( u) ]9 P
8. f1 = 49;
   6 u& p& i5 b, g) B
9. f2 = 51;
   : Z0 b% h* |$ ^! L. O9 i9 x
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;
    / x. A  q( _" x. k9 l
13. N = length(x);
    8 w* H: B& V1 y, D
14. pe = floor(N/2);
    
15. 
    
16. %% 构造样本矩阵
    # J2 g" {5 a! n& X5 ^' G
17. 
    . P; F8 x1 m8 A
18. Re=zeros(pe+1,pe+1);
    
19. 
    
20. for i = 2:pe+1
    ) F8 r! F% |/ X% Y8 l. o
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);
    
24. end
    & }$ u$ Q" R; r* h4 {
25. end
    ; a! d  _+ {6 k2 f% E
26. end
    
27. 
    ( D+ v# J) E7 g' c
28. Re(1,:) = [];
    1 ?' G( b- t7 a2 ^
29. 
    8 S4 l9 b, s1 u- m$ ]' O
30. %% SVD_TLS确定介数p及a
    
31. 
    4 m1 j, D" z" f) M' ^; w. M
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    6 q2 s9 a' {+ Y7 F  {6 @
34. % 求p值
    
35. % 计算全部奇异值平方和
    5 Z+ ?4 t* J2 ?
36. sum_all = 0;
    5 }0 U+ }# q( }& i
37. for i = 1:pe
    2 a4 }, x; W- w/ E7 C
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    
40. 
    , }* l. }4 P* ?6 G& i
41. % 归一化比值Ak/A 求p值
    
42. sum_k = 0;
    4 D4 d- V+ d* z3 _
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
    
48. p = k-1;
    % ]. |% P! `, [4 }/ ]) j; K
49. 
    
50. % 求Sp部分
    
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    
52. for j = 1:p
    % }( a, Z1 G4 v# n
53. for i = 1:(pe+1-p)
    
54. Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    9 Y3 E9 p/ z  G$ @! m
55. end
    3 d4 o. Q6 E, X9 X  m; ]0 B
56. end
    % d2 m/ b$ ~1 W# Y, }% u: f9 _. k
57. 
    
58. % SS = zeros(pe,pe+1);
    7 G" I. R% Y# G
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    ( c# Q( \) ?; E+ p% @
61. % end
    8 C& N& u# h7 d- l. F- f
62. % B = U*SS*V;
    
63. 
    8 @  P' `5 B: P$ z+ v
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
      S1 a7 L4 F' L7 g) N
66. if isinf(inv_Sp(1,1)) == 1
    3 U# z- {4 f" J- u6 A, L
67. inv_Sp = pinv(Sp);
    4 V7 }* i7 J7 E
68. end
    / h0 l  s9 d, F( X; c3 g9 j5 O
69. 
    & k" u' P4 [( V3 \
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    1 _4 B/ l% z0 m! E! c
72. 
    + m# J1 F2 J/ r" V: Z+ r6 U
73. %% 求z
    * I) Q- P& x- E5 c& [$ t
74. y=[1 a'];
    ! u3 B; {% L2 e2 F5 A
75. z=roots(y);
    
76. 
    
77. %% 求x的近似值x_j
    
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    
80. for i = 1:p
    + _& D8 p  O8 L; v- A: H
81. x_j(i)=x(i);
    7 `( d1 A, O6 N0 `; M. s! D
82. end
    2 _2 `& P! i" M* _' T3 U
83. 
    ; k9 A- t# R  n) h, _
84. %求x的N-p+1个近似值 x_j(p+1:N)
    
85. for n = p+1:N
    $ }% p% }2 F# V8 F6 p
86. for i = 1:p
    3 v: P" S1 K2 h
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    7 ^: W5 ^6 n1 o* \
88. end
    ) S; P* U* ^! _! S: L4 R
89. end
    
90. 
    
91. %% 画图 x、x_j
    ! J' ~: v% N2 M
92. hold on;
    ( ]7 T% e4 [) x5 L+ p  M
93. plot(t,x,'k');
    & w) z' n) }$ W/ @: W: r) _
94. plot(t,x_j,'r');
    8 ?: x4 X3 C/ u5 i8 u
95. hold off;
    , l  e# B) p6 l8 Y' b0 n
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    
98. 
    # w7 x' ]! X1 k  q+ Y# l+ l$ n
99. % 求取N X p维vandermode矩阵Z
    * X* o5 B% b' Z* G2 _0 {
100. Z=zeros(N,p);
     
101. for i=1:N
     
102. Z(i,:)=z'.^(i-1);
     2 V: m  D5 |- L7 g. T+ j/ p
103. end
     3 g, T3 Q6 V! S( d" V
104. 
     
105. %求取H
     
106. H=zeros(p,p);
     8 [- Q& b( E! i& ^  s  x3 ?
107. for i=1:p
     $ D- V( B; s+ o( i6 L/ Y3 L
108. for j=1:p
       b* |0 s) j8 D  @4 e: R/ @& F, {
109. m=(conj(z(i))*z(j));
       A: }" m! ]& l2 U# U4 O
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     
112. end
     ) }: R: }7 J$ H* W
113. 
     , r, W( W% g' L& A4 t
114. % 求取b
     
115. b=inv(H)*Z'*x';
     0 h$ V$ M' n% C: S# U
116. 
     - {# d/ f+ i9 a
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);
     ! C  K) o. s* P1 k/ g3 e; Y
121. Damp(i) = log(abs(z(i)))*dt;
     
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     / f0 Y! j0 y$ V/ {
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
+ a" l1 l4 p5 [( ]9 l2 X& f) Q

2 Q: l" c' T' e0 |* R    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006


! p9 H( T8 o2 \0 I- ^8 @    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
% ^+ [' ~1 o3 V0 t+ E, E( E/ a/ m
/ ?+ b* _! b/ Y
    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶