关于prony算法

xian2006

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

1. %% 数据准备
   
2. clear;
   
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
     e! L; Z3 U* ]4 x; w- \/ L
7. % t = t(201:500);
   
8. f1 = 49;
   
9. f2 = 51;
   " i/ O1 _$ K( E
10. t=0.0001:0.0001:0.01;
    9 _$ S# x* E2 h
11. x = 160*sin(2*pi*f1*t+pi/5)+150*exp(-3*t).*sin(2*pi*f2*t+pi/4);
    4 w/ c: b1 z4 [5 d, n2 @
12. dt = 0.0001;
    
13. N = length(x);
    1 ~3 l9 l2 \* p  J
14. pe = floor(N/2);
      w/ D' X# _9 i/ D! V( v% Q! H) a
15. 
    
16. %% 构造样本矩阵
    
17. 
    2 B- f" H' y8 D/ ]4 m+ z
18. Re=zeros(pe+1,pe+1);
    # e- h3 v4 X1 m1 q. a; |
19. 
    5 A; I* c7 U% G, `( r% @9 M; T
20. for i = 2:pe+1
    . N% @: r9 i# r# C3 B0 X$ o+ j
21. for j = 1:pe+1
    1 K/ B1 M* y( v/ G
22. for n = pe:N-1
    - G2 i( D1 o' ^: T7 s4 y: ~/ [$ _
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    ) L6 u  D+ ?  i* V
24. end
    / C+ @& q; d/ N1 {9 R8 \% U
25. end
    ' Z- B! I+ I  n) U. J+ U; X
26. end
    
27. 
      d. Q& \" R5 k* T  I7 z& R0 S
28. Re(1,:) = [];
    4 s/ X. b5 r% Q2 u& i. _
29. 
    
30. %% SVD_TLS确定介数p及a
    : [% o+ t# @8 z5 r& g
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    0 @9 \: q" o" S) L
33. 
    8 N; I7 o8 l. U) f' c0 |
34. % 求p值
    6 P# {! w' }3 {/ K" }) m: S
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    
37. for i = 1:pe
    + Y% l$ }9 {3 \9 [- @% z; e2 J
38. sum_all = sum_all+S(i,i)^2;
    
39. end
    0 ^! l& G; f7 {5 R0 m) f! c( l
40. 
    
41. % 归一化比值Ak/A 求p值
    + H/ j; r( }+ @' B- ^
42. sum_k = 0;
    
43. k = 1;
    
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    & M4 ]$ M! F4 [$ f% S
45. sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
    
46. k = k+1;
    
47. end
    ! h' U/ k) z% d6 n0 n0 k
48. p = k-1;
    
49. 
    0 J- D1 O7 D- C9 Y
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)
    
54. Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    
55. end
    
56. end
    ( H. q7 ~% f7 L8 h
57. 
    
58. % SS = zeros(pe,pe+1);
    
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    & e* M9 |* R3 }3 T/ z+ \6 v
61. % end
    3 n+ z: ~, G$ \. N
62. % B = U*SS*V;
    
63. 
    $ n( u: z, O/ O: _0 b6 O0 Q: _! O
64. % 求Sp逆矩阵
    
65. inv_Sp=inv(Sp);
    " @! K3 z8 n  w4 I
66. if isinf(inv_Sp(1,1)) == 1
    
67. inv_Sp = pinv(Sp);
    ( o' m- e/ s% q% _  o; S( b" h
68. end
    " Q6 }/ f- x. s, w3 |) e; L; Z
69. 
    
70. % 求a
    6 w8 J3 D- L3 Z! g! C3 T  U0 g
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    
72. 
    1 P# k) ?  O8 D5 _
73. %% 求z
    
74. y=[1 a'];
    
75. z=roots(y);
    
76. 
    % K# ^! q7 z, {  j6 m2 i" Y9 `9 O
77. %% 求x的近似值x_j
    % ~8 |: L3 d' r% _( }! O2 ^# z% R
78. %求前p近似值等于测量值 x_j(1:p)
    & U& {  I- l% f& Z% {2 v" i" |
79. x_j=zeros(N,1);
    
80. for i = 1:p
    
81. x_j(i)=x(i);
    
82. end
    6 w# d* Q6 x5 @9 v% G0 {3 }
83. 
    8 O; r5 [3 q3 E/ C" K# Y
84. %求x的N-p+1个近似值 x_j(p+1:N)
    & S7 k. U% l( Y+ K, q
85. for n = p+1:N
    / v6 S, H) U5 K/ j9 y: f
86. for i = 1:p
    
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    $ g4 D4 e5 u3 k9 {0 O0 t2 A
88. end
    
89. end
    # y. U1 h: @% G: D4 q8 J& `( \
90. 
    5 C2 t! U! o/ O) e
91. %% 画图 x、x_j
    " C+ A; C( h* `( {7 n
92. hold on;
    
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    
95. hold off;
    9 i2 H* M5 f9 c& Y
96. 
    
97. %% 求取 b=inv(H)*Z'*x_j
    8 _$ \  P2 c2 Y/ B$ [
98. 
    $ D: R: `- {" x6 {9 B% i% ^
99. % 求取N X p维vandermode矩阵Z
    0 O2 p- M& M& D' e
100. Z=zeros(N,p);
     ) o+ w. H  Y; R. f' N7 Z
101. for i=1:N
     
102. Z(i,:)=z'.^(i-1);
     ! d3 s( R3 b2 m' J5 X; Z8 f/ F
103. end
     
104. 
     
105. %求取H
     
106. H=zeros(p,p);
       ?) ^  U  n3 o8 z
107. for i=1:p
     
108. for j=1:p
     6 D) Z: p! |+ i7 n8 r( O5 X% m
109. m=(conj(z(i))*z(j));
     6 t3 T# b7 `; F/ N
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     1 I5 U( h+ `( T2 K! J! k
112. end
     
113. 
     
114. % 求取b
     1 D. C- L% F) ]0 \3 z( g! p
115. b=inv(H)*Z'*x';
     
116. 
     9 m. `( \7 \+ J( V' V& `
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     5 Q4 h: }, k* q1 i& n/ f
119. Amp(i) = abs(b(i));
     5 M; W9 s1 y' k
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))));
     
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
7 S$ l6 B7 U# k$ c  |

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

sjx0323

回复 7# xian2006
# v$ \( q# T. @  E

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

sjx0323

回复 7# xian2006
0 l( n; w/ S! f4 V- G

) I$ K* I" v: t$ ]8 f4 e0 t/ q    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶