关于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的原程序，发现工具包里边的算法和上边两本书上的算法是不一样的，好像是根据零极点和滤波脉冲响应什么的，我实在不是太懂。
/ P" b. p! x0 M
  j) T9 D9 W3 {  D所以在这里诚心向各位请教

1. %% 数据准备
   
2. clear;
   
3. clc;
   ( l1 O2 g* S9 z# f7 \& r5 m
4. format long
   , v0 X' I2 @, ^) F7 I
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   7 }' P. g8 \4 b2 [  B0 Z8 D: Z6 e3 }
7. % t = t(201:500);
   " |5 a. f. q+ W% R% b
8. f1 = 49;
   7 }- }4 `- r9 m7 J/ b  L6 E( Y
9. f2 = 51;
   * k6 s8 Z1 W) Z, U, L+ C' \1 a
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);
    % N/ B/ n9 @; z. w
12. dt = 0.0001;
    1 M4 N$ A0 o! ~, ?4 R
13. N = length(x);
    
14. pe = floor(N/2);
    9 _; I* X5 S; b
15. 
    ( J8 Z4 B! a4 M; O$ v5 M8 j" Z
16. %% 构造样本矩阵
    
17. 
    8 e$ r" v% _" g. P2 y. ?
18. Re=zeros(pe+1,pe+1);
    , y( i( j* r9 j7 W# c, W
19. 
    
20. for i = 2:pe+1
    ' B0 C% o) N7 f1 ^
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);
    , s5 l" a9 r' |
24. end
    ' S& X1 ^% l/ y0 D, x  {
25. end
    
26. end
    
27. 
    - E- k; U* K! _" l3 R0 d' k
28. Re(1,:) = [];
    
29. 
    ! ?: J/ |4 s$ V( M6 _
30. %% SVD_TLS确定介数p及a
    3 B9 r2 {3 c/ k/ X$ R1 d3 ^
31. 
    8 Y" |4 Z1 i3 o: D, P. r! `
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    8 c5 Z3 G; E* X5 m2 S
33. 
    + ]/ p3 y$ V7 q
34. % 求p值
    . E! v8 h: b2 M% w: m4 I& ?0 V2 V+ b
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    1 o* W+ P( k) R& W, z) X
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    
40. 
    + _' R  n1 i) r8 U
41. % 归一化比值Ak/A 求p值
    
42. sum_k = 0;
    
43. k = 1;
    2 o3 @$ d1 N( R, t
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;
    - H  J) j/ U, {+ [  }% X
47. end
    : z; t8 _; y8 U# R/ h- t7 q" N8 o
48. p = k-1;
    & g9 o$ s1 T, ~! P% [( D* i1 @
49. 
    + }6 e8 M' ~; l1 F' l( C- b
50. % 求Sp部分
    
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)
    $ W  M# K9 h1 V, a3 e  |( T
54. Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    1 c4 ~( B2 T' g3 m" r
55. end
    3 N1 D( U8 Z/ ~+ @
56. end
    ; e5 E; q6 _. O* q- m6 C5 J# I7 o
57. 
    
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    
61. % end
    / q% g7 ^- q2 E, ?) m9 M% G
62. % B = U*SS*V;
    
63. 
    7 Y& q4 d, p/ M
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
    ) ]+ m* Z6 @$ o( ?9 Q% }
66. if isinf(inv_Sp(1,1)) == 1
    
67. inv_Sp = pinv(Sp);
    $ `3 j6 @5 n* o! @+ r7 G
68. end
    
69. 
    
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    ( k  y, X1 H( @1 V1 z/ H
72. 
    - Z9 j3 s9 q7 w. g6 i- s9 O
73. %% 求z
    
74. y=[1 a'];
    
75. z=roots(y);
    , \! p+ c9 l4 l- c: D% J9 t: z
76. 
    . M1 Q$ ~8 V& Z4 d+ _! U6 @$ l) C
77. %% 求x的近似值x_j
    ! g9 Q& g6 X  S: M
78. %求前p近似值等于测量值 x_j(1:p)
    2 @$ A. s) M: t- f8 j, Y
79. x_j=zeros(N,1);
    / u5 I4 X& N# \4 d* E4 v
80. for i = 1:p
    
81. x_j(i)=x(i);
    " b2 F' I( ]7 J& h8 d$ S. h
82. end
    
83. 
    
84. %求x的N-p+1个近似值 x_j(p+1:N)
    ' X9 ~! S5 M" k* `
85. for n = p+1:N
    
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    ' G6 x$ b/ R, h6 \/ w& S# m
88. end
    
89. end
    6 \7 q( c% t- z. A
90. 
    5 U3 C) q, o9 t
91. %% 画图 x、x_j
    
92. hold on;
    7 _: c6 }) o8 t
93. plot(t,x,'k');
    ' K5 I* l! e4 a+ H
94. plot(t,x_j,'r');
    
95. hold off;
    
96. 
    6 R1 M7 ^# b# X& |
97. %% 求取 b=inv(H)*Z'*x_j
    7 T0 w1 E/ z- I- u
98. 
    9 `) ?5 O' X2 P* o
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     
101. for i=1:N
     ( r1 L" D3 G. c% a
102. Z(i,:)=z'.^(i-1);
     % F9 @4 i6 {. {
103. end
     & T+ K! X9 v1 D0 @0 O
104. 
     
105. %求取H
     
106. H=zeros(p,p);
     8 j" R9 _  Z& H: Z. D6 M
107. for i=1:p
     ; \; k/ |" N5 g+ @9 F$ ]& L
108. for j=1:p
     9 X) \. u# \2 z2 h) R
109. m=(conj(z(i))*z(j));
     
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     
112. end
     ! t2 r: j1 |& A/ d
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     7 Q! H6 A8 j+ _: B0 f
116. 
     1 ]5 t, z0 {' V
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     
119. Amp(i) = abs(b(i));
     ! O/ ~! ~) |" @( B; m! v
120. Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     
121. Damp(i) = log(abs(z(i)))*dt;
     ) R% g2 c9 v) [7 r+ j
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
+ v) ]( h: f; L) k) S) j, n

    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
; c3 G1 [; t  `0 G/ F7 K

, K9 ^% I* }: L5 h& @5 o/ B: z    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
: ~2 v) _. A) v( |8 L: Y" J9 B
/ l% y9 A( X% A5 q. E
    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶