关于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的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
3 Q4 U; A+ Z( W- ?) m
0 o4 Q% G( N* {所以在这里诚心向各位请教
/ g+ k3 ~& I+ d% r

1. %% 数据准备
   % S3 p" U/ j. p
2. clear;
   # z, p4 N, r- g# y( ~! Z% c3 E
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   ' j! G! f& {$ ]# z% u8 D6 X- O
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   0 y( D. N6 o2 n$ e+ g
8. f1 = 49;
   
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);
    4 B  |  n) w( f  q+ F9 u
12. dt = 0.0001;
    $ {7 Z) v& r' d0 C# v, o
13. N = length(x);
    
14. pe = floor(N/2);
    * b9 l  ^0 ?7 j
15. 
    
16. %% 构造样本矩阵
    2 S5 t+ \( _% X; ~3 o3 ?0 A
17. 
    
18. Re=zeros(pe+1,pe+1);
    : x% L. d" a9 T* i0 I
19. 
    
20. for i = 2:pe+1
    2 j4 r& i: v" b8 z2 Q
21. for j = 1:pe+1
    : Y3 s# ?$ ?) P6 I6 {- M4 [
22. for n = pe:N-1
    
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    
24. end
    , W0 O- _! ~. |
25. end
    1 m9 t4 Y; y( D: ]/ Q  N8 o3 ]# r
26. end
    ' F* F8 s$ G. z5 z
27. 
    / t; w* b6 B0 P" Y, i
28. Re(1,:) = [];
    0 R+ u3 ?+ L9 n( G+ o
29. 
    
30. %% SVD_TLS确定介数p及a
    
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    : P' q6 O3 _$ }6 }* g0 |
33. 
    1 U* O: a9 a/ V/ X  q4 q
34. % 求p值
    " N) z, R# T* C# t
35. % 计算全部奇异值平方和
    : v+ \6 Q8 P5 }/ u
36. sum_all = 0;
    $ x0 F) a, V7 C) B, o
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    
40. 
    
41. % 归一化比值Ak/A 求p值
    
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个奇异值平方和
    - G& L& x! v! d7 l  S) D' ^; S( ^4 j$ T4 A1 C
46. k = k+1;
    
47. end
    6 H$ l0 g% v3 d1 e& L( o; K) C: g' i
48. p = k-1;
    
49. 
    $ O6 H- L; X3 M% C
50. % 求Sp部分
    
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    7 O6 f8 C" ?2 @0 U' Z; c6 O
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
    , i" X! l# t1 z% O( X
56. end
    
57. 
    
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    8 t, K; x- o& `7 V0 _. G
61. % end
    * a! ?- w8 F0 x8 ]3 y
62. % B = U*SS*V;
    % }. R" n+ _3 ?. l3 r
63. 
    : O6 L; p' \0 X3 u5 H: r
64. % 求Sp逆矩阵
    7 S: B3 H2 h4 o& x6 h- ?
65. inv_Sp=inv(Sp);
    . h- m% @6 A- T4 v0 c
66. if isinf(inv_Sp(1,1)) == 1
    6 h( i: h! n. S6 m& s2 B* O# y
67. inv_Sp = pinv(Sp);
    , Q. N! b* [2 }
68. end
    
69. 
    3 A" L. B2 `* ^, k2 a
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    % Y% T+ @+ p' M( A; }6 R
72. 
    
73. %% 求z
    , a) \$ H; `& w# p/ {
74. y=[1 a'];
    
75. z=roots(y);
    
76. 
    . C- k5 ]+ K( s
77. %% 求x的近似值x_j
    2 f* X0 _2 u3 B# j6 l  Y3 D2 ~
78. %求前p近似值等于测量值 x_j(1:p)
    3 l7 }) K& @0 Q& J$ Y6 F+ I
79. x_j=zeros(N,1);
    
80. for i = 1:p
    4 _) z/ X% J. U8 L, f+ w# O! k
81. x_j(i)=x(i);
    % L9 Z& D2 A/ s( T# Y
82. end
    
83. 
    ' e5 s( j1 a8 _' S4 i( X
84. %求x的N-p+1个近似值 x_j(p+1:N)
    ! S* }" o. W' l7 I- b! x( w
85. for n = p+1:N
    
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    , E2 z8 F; r( d! C. }; }4 c  ?  U
88. end
    8 b( v6 s5 M/ m5 {# R3 H; q; j
89. end
    , l" ]* I( W! H+ D
90. 
    
91. %% 画图 x、x_j
    " t5 H9 l4 j8 s" _
92. hold on;
    " N/ o: p4 ^: w$ s/ ^3 ?/ H
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    , Y% i* x6 q% [& m  M
95. hold off;
    
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    
98. 
    
99. % 求取N X p维vandermode矩阵Z
    : i4 l8 Z) Z) ^
100. Z=zeros(N,p);
     1 d1 B  `0 }6 q
101. for i=1:N
     0 X& D* I5 o6 R: R2 y. `
102. Z(i,:)=z'.^(i-1);
     " c; a, A3 E/ v! R
103. end
     
104. 
     ( N; }7 [+ r  J) J% _$ a4 J
105. %求取H
     . i) g+ I/ Y% i( W
106. H=zeros(p,p);
     
107. for i=1:p
     
108. for j=1:p
     
109. m=(conj(z(i))*z(j));
     9 q2 Z" F3 A% ?% p9 i1 Q
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     # O. P% s$ M5 q% B$ n
112. end
     4 o6 M' A; \0 K1 D) k9 W
113. 
     , O, K- y7 ^: {' L: m% a; P' F
114. % 求取b
     
115. b=inv(H)*Z'*x';
     
116. 
     : I* @& g: z3 b4 F
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     2 q, Q& T. P' o2 |. e, Q" H. C" r
119. Amp(i) = abs(b(i));
     6 p1 t/ {) t' t. r: x: c
120. Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     9 i. H# K/ E# Y* a
121. Damp(i) = log(abs(z(i)))*dt;
     6 |4 F7 @9 s. A* }/ K, q# }0 P
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     $ P" h, I5 I! {! e1 b6 f
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
4 n9 C2 f, l6 o( u+ m. u
2 A5 e1 a. p, ~  w; V; J
    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
0 O$ D) n$ F6 M

. K" z8 |1 {: q; h7 f5 m- y    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006

- A/ j  z  s/ L: S0 j
1 k6 F  f* ?" I    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶