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

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

1. %% 数据准备
   ) G% A* ^' J( b8 b" o% ~# o
2. clear;
   ! U- _3 {8 h7 _" s7 A7 _- C
3. clc;
   9 o, H, Y3 r2 U( @" M! u6 g1 |* k: D
4. format long
   2 l0 D( b! q: P- ^
5. % load('1000kV示范工程线路','t','vX0043a')
   4 y  M7 B/ a0 e
6. % x = vX0043a(201:500);
   & s$ D; D7 F5 [1 l! h9 u6 |/ s
7. % t = t(201:500);
   7 r" _- o' E$ b
8. f1 = 49;
   $ ]* H! y( ]3 \5 |  {
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);
    + Q6 f# Y1 R& C. f
12. dt = 0.0001;
    4 u, @# O% \  A! H9 b, {
13. N = length(x);
    ) o8 e- R! R- X8 e% O
14. pe = floor(N/2);
    $ L5 ~+ @! p' Y: Y. c" H
15. 
    2 m. d! G8 _# H
16. %% 构造样本矩阵
    
17. 
    % T- g) w8 q0 m0 x  l$ I" k3 h/ q
18. Re=zeros(pe+1,pe+1);
    
19. 
    : k9 Q+ i3 x1 b( k9 P+ ?. ~
20. for i = 2:pe+1
    
21. for j = 1:pe+1
    * x" F  ^& n/ F: Y
22. for n = pe:N-1
    1 q: Y% m: }6 n; J; q- _! ~9 X
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    % |# n' `3 q, X3 Q6 Z! x
24. end
      V2 {/ G! M3 E( A5 }
25. end
    5 _- U: ^1 v8 S3 k, w) X$ U0 ~
26. end
    
27. 
    1 C, n1 H$ T7 |* R4 i" m3 j
28. Re(1,:) = [];
    * ?, t: K$ l9 e9 A2 x# E4 L
29. 
    7 ]' n$ l5 Y! Q3 [# P5 @& q- @
30. %% SVD_TLS确定介数p及a
    5 \* E( Q9 e- [: N
31. 
    
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    
34. % 求p值
    
35. % 计算全部奇异值平方和
    
36. sum_all = 0;
    
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    ' C2 z6 T* n. D$ D
39. end
    
40. 
    
41. % 归一化比值Ak/A 求p值
    , C6 J+ U$ `( ~
42. sum_k = 0;
    . ?: v$ ^8 O7 _8 p
43. k = 1;
    
44. while 1 - sqrt(sum_k/sum_all) > 0.0000000000000000001 & i<=pe
    ( Y  U( f- R& s$ ~( @3 r
45. sum_k = sum_k+S(k,k)^2; %%%%%%%计算k个奇异值平方和
    . {! D- E! S/ n; T" h
46. k = k+1;
    3 k- b; W6 q$ [6 c
47. end
    
48. p = k-1;
    : k, [! F& M; m8 x3 e- {
49. 
    0 `4 k! Z% w$ g# D# C: {* r  a2 O* A
50. % 求Sp部分
    
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    8 x7 J" ~  ^. s( q
52. for j = 1:p
    5 w8 E+ o0 `( r. x5 x) g# h
53. for i = 1:(pe+1-p)
    
54. Sp = Sp+S(j)^2*V(i:i+p,j)*V(i:i+p,j)';
    4 }5 H; U$ e, e, [
55. end
    
56. end
    " z* o7 H7 |! ~" T
57. 
    & Z1 w" E) p8 g5 W+ F( c- [# m
58. % SS = zeros(pe,pe+1);
    5 o. z! |- F* D+ T1 V
59. % for i = 1:p
    5 m3 W" A7 [  a8 y+ k% H
60. % SS(i,i) = S(i,i);
    
61. % end
    % W# K9 ~0 z9 P# x7 H- a
62. % B = U*SS*V;
    ; c  r4 F  ?& G4 ~. R9 V8 l
63. 
    
64. % 求Sp逆矩阵
    % d! \) l8 e4 o- M; {8 x
65. inv_Sp=inv(Sp);
    
66. if isinf(inv_Sp(1,1)) == 1
    
67. inv_Sp = pinv(Sp);
    1 L+ z, K' R9 j" _) ], O
68. end
    
69. 
    
70. % 求a
    
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    
72. 
    
73. %% 求z
    ) i0 l( _) j( m9 N
74. y=[1 a'];
    ! e& S/ U5 ^# C6 Q' u  o6 |/ @
75. z=roots(y);
    ; J/ n: o# J- I$ G) x3 f& V
76. 
    $ x' T. \/ [) v+ L  x, J$ Z' y$ C
77. %% 求x的近似值x_j
    
78. %求前p近似值等于测量值 x_j(1:p)
    
79. x_j=zeros(N,1);
    
80. for i = 1:p
    
81. x_j(i)=x(i);
    8 Y: _0 f' F: V4 O% O
82. end
    : K- O5 Y0 |* t- H) A: w
83. 
    
84. %求x的N-p+1个近似值 x_j(p+1:N)
    0 _# @' f0 ~4 |- {! ]" i
85. for n = p+1:N
    
86. for i = 1:p
    . j0 X; F8 [" p4 a' C' |1 P1 F
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    
88. end
    ) ]& f, {7 _2 o8 |/ s7 Q
89. end
    
90. 
    / F" i% D6 q. b) a1 p. M
91. %% 画图 x、x_j
    - s' E* f4 G- ~8 W6 h+ f, ?
92. hold on;
    * E  b- I, o3 G+ M
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    
95. hold off;
    
96. 
    & P7 K" t( [- h9 Q  g: u
97. %% 求取 b=inv(H)*Z'*x_j
    + x9 F: z- R3 `4 R- p
98. 
    
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     
101. for i=1:N
       ^, [7 u. V7 _8 `" {
102. Z(i,:)=z'.^(i-1);
     * O. O% ]; s; a6 G7 y$ h2 z8 R
103. end
     $ M: |, a9 [7 Z) j
104. 
     # G6 A7 v* L( `. x
105. %求取H
     + e8 [  ~# F) v& d1 }1 y
106. H=zeros(p,p);
     
107. for i=1:p
     
108. for j=1:p
     
109. m=(conj(z(i))*z(j));
     + ^3 K. j! A: h7 T+ y
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     
112. end
     
113. 
     
114. % 求取b
     
115. b=inv(H)*Z'*x';
     
116. 
     
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     
118. for i = 1:p
     6 G; }- ^  R. G# T' H" ]
119. Amp(i) = abs(b(i));
     / u. n+ H8 ~8 t; ^4 Q
120. Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     
121. Damp(i) = log(abs(z(i)))*dt;
     * x+ i3 f- z% G
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     4 H" \; u  j3 ^7 ?6 \
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323


' _4 Z/ S4 @0 d6 B3 [. k    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
  c) o* p( O, o- V$ n

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

sjx0323

回复 7# xian2006

; _% R; V5 r  A8 V1 M2 d5 _  R2 q
    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶