关于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的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。

所以在这里诚心向各位请教
# S" Y) J- z' k' A" `

1. %% 数据准备
   
2. clear;
   . m/ ^) f& ]$ N6 D3 N& e
3. clc;
   9 @/ U9 |* ?8 @; U" ?8 e; ?
4. format long
   / K# ?0 C9 ^' ]. T2 c& \* `4 T
5. % load('1000kV示范工程线路','t','vX0043a')
   9 @7 f2 ^& z( j" x) I% E
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   ! X) R  ~' w. k( i- p! T- s
8. f1 = 49;
   . u) M, q  d+ v5 Q
9. f2 = 51;
   2 K6 [1 b# x$ g0 ?
10. t=0.0001:0.0001:0.01;
      W9 {: E4 Q  Q) X* G& X! T
11. x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
    " Y9 T2 t$ p4 j; m3 n
12. dt = 0.0001;
    2 r$ m( _- {+ d6 j. j/ b
13. N = length(x);
    
14. pe = floor(N/2);
    6 n9 W: Y9 x0 P- J- ?
15. 
    1 i, @) T3 |4 {# K( v
16. %% 构造样本矩阵
    7 E  l3 u9 i4 D# b6 \
17. 
    
18. Re=zeros(pe+1,pe+1);
    
19. 
    
20. for i = 2:pe+1
    * [$ A* t1 j( t( ?
21. for j = 1:pe+1
    # T# U# f6 s$ M8 u" o
22. for n = pe:N-1
    
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    
24. end
    , _, U6 a1 r6 B7 W* X$ \& ^
25. end
    - E" P5 \: o' B7 f, r
26. end
    % I; T5 k+ D: }8 D3 Z4 Y# b
27. 
    
28. Re(1,:) = [];
    
29. 
    $ b  m; z5 |, @+ K' e) Q
30. %% SVD_TLS确定介数p及a
    5 N: Y" {/ s' G; H$ ], v4 @7 {
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    * B# Y6 F/ i0 ~
33. 
    
34. % 求p值
    ) A1 W. n5 t  c% L% p$ q
35. % 计算全部奇异值平方和
    ( A! ]7 E+ X9 a1 z
36. sum_all = 0;
    
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    ; {- e0 Z! n) P/ u! Z
40. 
    0 t$ [. K/ B# P3 X2 z6 ^2 L, U  H$ S
41. % 归一化比值Ak/A 求p值
    ! m* N9 W7 J. R  \' @5 g. e
42. sum_k = 0;
    6 W/ E9 F5 P1 ]; o
43. k = 1;
    1 @+ v1 m+ M- `* l' |; r/ R8 d( `
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    ' U2 g* f$ m0 u! ~8 K6 q) I
45. sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
    
46. k = k+1;
    5 v! R0 D( Q. V8 @  V' r! W; E
47. end
    2 ?, B1 Q2 F) |7 t7 j+ w
48. p = k-1;
    
49. 
    
50. % 求Sp部分
    ) T! |$ w6 o/ G5 X4 Z
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    
52. for j = 1:p
    9 k* r4 r' T6 _1 ?  h3 i
53. for i = 1:(pe+1-p)
    ; v: H2 N/ k! Z! H% }4 h
54. Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    # J/ M- i8 U4 _; k, S3 z
55. end
    $ U# X1 P* B* t) L
56. end
    9 ^$ T9 j7 t4 q( |$ K/ ]
57. 
    $ L" t9 J0 p' H( l
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. 
    " ?* n: B/ S$ z3 o7 F
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
    
66. if isinf(inv_Sp(1,1)) == 1
    ' u2 Y5 b3 R2 J& w
67. inv_Sp = pinv(Sp);
    9 g3 y+ X# R) Q  w( a: t/ v5 a, _
68. end
    
69. 
    
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    * _& Y0 a: p# P8 ?: ~6 i  ~
72. 
    7 F1 L" \7 T; ^' G
73. %% 求z
    1 G3 o" Z0 o; B
74. y=[1 a'];
    
75. z=roots(y);
      L3 `- n  M! W
76. 
    
77. %% 求x的近似值x_j
    & h1 u; @& O  I1 W
78. %求前p近似值等于测量值 x_j(1:p)
    ! z. `  y1 |3 i7 a
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)
    * y+ X+ s  J4 N: Q& ^
85. for n = p+1:N
    
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    
88. end
    ' \* z. H# _: H/ V' d
89. end
    : ?+ e" f6 m0 d% T2 d
90. 
    8 K* R" s" |" h/ E9 ^
91. %% 画图 x、x_j
    2 p, u9 ]# s* C0 @
92. hold on;
    
93. plot(t,x,'k');
    ) U; }: n9 B: V- j2 j" ~1 o& i
94. plot(t,x_j,'r');
    / S/ q/ R0 {3 Q' N' _& l
95. hold off;
    ' _+ R) ]8 j2 A4 b* s
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    
98. 
    
99. % 求取N X p维vandermode矩阵Z
    + \- E& j3 H6 `+ Z$ U
100. Z=zeros(N,p);
     , P; K: c- T* c& d
101. for i=1:N
     1 p9 r& V* T: G
102. Z(i,:)=z'.^(i-1);
     7 b) x+ O3 ^: a* H6 y; h) f$ r9 d+ V
103. end
     
104. 
     
105. %求取H
     
106. H=zeros(p,p);
     7 J3 _2 K9 q( R- a- P
107. for i=1:p
     3 i# @; _9 G7 L
108. for j=1:p
     ) A4 }' e2 j5 G/ \: j
109. m=(conj(z(i))*z(j));
     
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     % k% }# `2 h8 D; ~) X
112. end
     , y2 c/ i! k% x/ Y- j& f  g4 X- H
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     $ \) I* l+ J& ^& H
116. 
     6 i$ g" T( y, P7 r* J, ^
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     ) z2 D) q2 q) U$ D
118. for i = 1:p
     4 k+ O/ r6 M' z4 O9 ^
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;
     
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     ' p, h( \1 s4 b" Y+ a
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
" ]7 g: B" \7 J; v0 \9 ~
: D1 m4 z- d; O9 H0 p; p
3 L6 }- T  y5 k, `    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
7 Q* v/ x9 e6 v/ ]' S' J

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

sjx0323

回复 7# xian2006
+ V5 s6 b) B( [. n
7 n% W% b2 T$ W% T$ {( O
! S9 R# L  h# O* v# i    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶