关于prony算法

xian2006

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

1. %% 数据准备
   
2. clear;
   
3. clc;
   
4. format long
   
5. % load('1000kV示范工程线路','t','vX0043a')
   
6. % x = vX0043a(201:500);
   
7. % t = t(201:500);
   ' d, A" L# K2 I4 }
8. f1 = 49;
   7 g$ @3 g# z( X* z' }
9. f2 = 51;
   ) ~+ o9 |3 E2 j. d
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);
    9 c& S2 F3 q! _% K
12. dt = 0.0001;
    4 A: [" s: O" V
13. N = length(x);
    ! U* a; a. O8 i
14. pe = floor(N/2);
    8 W. f, |; `+ d7 v4 n$ Y
15. 
    2 E3 h  B- G5 p/ L" r, H
16. %% 构造样本矩阵
    : G* V% B. p2 I
17. 
    
18. Re=zeros(pe+1,pe+1);
    
19. 
    & u3 U+ J2 [  n
20. for i = 2:pe+1
    0 p+ P2 A2 r/ Z  K& j
21. for j = 1:pe+1
    
22. for n = pe:N-1
    3 m7 b' V. ^) n- i( z1 S) w/ ^7 n5 ?+ @
23. Re(i,j) = Re(i,j)+x(n-j+2)*x(n-i+2);
    
24. end
    # Z/ f) q4 R, v5 d- s
25. end
    
26. end
    9 s- R, ~% q' u7 A0 n) L' M
27. 
    
28. Re(1,:) = [];
    
29. 
    
30. %% SVD_TLS确定介数p及a
    9 v' D; P' P( f
31. 
    % G, Z* s( X* A
32. [U,S,V] = svd(Re); %%%%%%奇异值分解
    
33. 
    
34. % 求p值
    1 n1 H- k7 \* y( i5 X5 |- V
35. % 计算全部奇异值平方和
    3 k4 I2 \* ^) M6 G3 _5 ]
36. sum_all = 0;
    6 O- ]' s% V3 j/ Y* L$ A
37. for i = 1:pe
    
38. sum_all = sum_all+S(i,i)^2;
    * E& N5 b0 i( ?
39. end
    
40. 
    
41. % 归一化比值Ak/A 求p值
    2 Y+ |" z) r) ^- u1 f: B; c+ @7 k; I) C
42. sum_k = 0;
    
43. k = 1;
    ! M( z3 R0 @, a1 K& P1 P
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;
    8 S+ q% B: ?6 o* |& Z
47. end
    
48. p = k-1;
    - R& o; W- R9 N5 M
49. 
      o& D  F) P5 s) x, t
50. % 求Sp部分
    4 h+ Y9 K# t9 n4 q' A! t
51. Sp=zeros(p+1,p+1); %%%%%%s生成（p+1）X（p+1）维矩阵Sp
    * p% t7 k6 b  ^
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)';
    ( g/ J2 S3 z* \7 d3 t/ m. |
55. end
    
56. end
    . o4 W( u: F4 q' A& s4 L9 {
57. 
    
58. % SS = zeros(pe,pe+1);
    # A4 O- |) Q9 m+ ^
59. % for i = 1:p
    
60. % SS(i,i) = S(i,i);
    
61. % end
    * s% Q( r9 ~* p
62. % B = U*SS*V;
    
63. 
    
64. % 求Sp逆矩阵
    ) _9 P$ {: @, ?* g0 M% r
65. inv_Sp=inv(Sp);
    
66. if isinf(inv_Sp(1,1)) == 1
    1 H; c8 T/ h" E6 Q& l, |" C
67. inv_Sp = pinv(Sp);
    + ]( [4 y0 k/ H4 V" \3 g
68. end
    
69. 
    
70. % 求a
    5 ^# L9 H' q- j- w
71. a=inv_Sp(2:p+1,1)/inv_Sp(1,1);
    9 X" s' G' p2 v; g1 S4 S
72. 
    2 D% N! V6 H9 `# m! S: y0 \
73. %% 求z
    
74. y=[1 a'];
    
75. z=roots(y);
    
76. 
    
77. %% 求x的近似值x_j
    
78. %求前p近似值等于测量值 x_j(1:p)
    0 {: f! D4 |% V2 X, T: ^: u
79. x_j=zeros(N,1);
    
80. for i = 1:p
    
81. x_j(i)=x(i);
    
82. end
    % E; q0 O. w( b$ Y
83. 
    3 z8 x& ]7 P. ?5 J
84. %求x的N-p+1个近似值 x_j(p+1:N)
    
85. for n = p+1:N
    
86. for i = 1:p
    / l* s5 z; q' H3 ?# t
87. x_j(n)=x_j(n)-a(i)*x_j(n-i);
    
88. end
    5 a! X. ^" e' S
89. end
    
90. 
    1 O9 q+ J" R3 [. D: l
91. %% 画图 x、x_j
    0 n1 Q; g# e) ]9 D$ I
92. hold on;
    ; r- g8 ^, j, U' D4 Y1 J
93. plot(t,x,'k');
    
94. plot(t,x_j,'r');
    ) @' V$ S- w) V3 h0 W; g
95. hold off;
    
96. 
    / p" s2 Z; V7 o" I8 O: U4 ~
97. %% 求取 b=inv(H)*Z'*x_j
    + B% t( e. F6 n6 V4 |$ C1 R- x
98. 
    
99. % 求取N X p维vandermode矩阵Z
    
100. Z=zeros(N,p);
     
101. for i=1:N
     - ~# J) ]1 s7 G
102. Z(i,:)=z'.^(i-1);
     
103. end
     8 U' F$ g& B# C( I7 b% x3 F
104. 
     
105. %求取H
     3 ~, c$ A) c! B) e1 q  o2 C& l+ w" [
106. H=zeros(p,p);
     
107. for i=1:p
     
108. for j=1:p
     " G  |3 m# F; N* k' |
109. m=(conj(z(i))*z(j));
     1 P/ V! Z' e0 f: x9 i
110. H(i,j)=(m^N-1)/(m-1);
     
111. end
     ! E; C# @. T+ y
112. end
     ; }& p5 i7 `8 Z" m
113. 
     0 E* O3 M2 {" B6 M8 R5 m
114. % 求取b
     3 y* a8 K0 {- E5 J$ ~+ {6 U# Q1 S3 ^
115. b=inv(H)*Z'*x';
     ! h0 `+ |8 q: t9 o8 [$ o
116. 
     ! i0 O0 q4 X6 t5 @; S
117. %% 计算振幅Amp 频率Fre 衰减因子Damp 相位 Pha
     3 g2 \) E* \3 }' Z( Q3 ^% _
118. for i = 1:p
     ! y/ @5 G" q8 ]% F/ [
119. Amp(i) = abs(b(i));
     1 L- Y2 B( f2 X3 s9 e* S5 b$ X& @
120. Fre(i) = atan(imag(z(i)/real(z(i))))/(2*pi*dt);
     0 p( E$ ^8 [) I$ C& M0 B
121. Damp(i) = log(abs(z(i)))*dt;
     9 Q  J$ @! H3 M; r) h
122. Pha(i) = atan(imag(b(i)/real(b(i))));
     ) c& k4 j9 f* A; g$ Q
123. end

复制代码

李真

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

Deallee

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

xian2006

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

tianlangistl

好专业，看不懂帮顶

sjx0323

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

sjx0323

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

xian2006

回复 6# sjx0323
  B, e: p% c6 |) q4 R$ X
% K; D/ P( q+ }
    真的吗 那真是太谢谢了 可以联系我吗？qq345749437

sjx0323

回复 7# xian2006
: T, L2 r4 e; H1 p/ Y: S* H
6 x4 |) n$ G7 o- x. k+ i
8 @2 q6 Q- Z. s5 V, p8 w    他今天不在，明天来了帮你请教一下！

sjx0323

回复 7# xian2006
& l* z9 k* b7 F# ~9 @) P

    看来这个问题解决不了，我似乎发现你貌似就是我的学长啊！呵呵

nanshugong

我帮学长顶顶