用Burg算法求AR模型的参数

VibInfo

主程序：

1. c----------------------------------------------------------------------
   
2. c   main program harburg : to test subroutine ARBURG
   
3. c   To compute the AR model coefficients by Burg algorithm ,IP=10
   
4. c   Please  link subroutine ARBURG,AR1PSD RELFFT,PSPLOT
   
5. c----------------------------------------------------------------------
   
6.         complex x(0:127),a(0:10),ef(0:127),eb(0:127)
   
7.         data n/128/,ip/10/,mfre/4096/,ts/1./
   
8.         open(3,file='test.dat',status='old')
   
9.         do 20 k=0,n-1
   
10.            read(3,*)x(k)
    
11. 20      continue
    
12.         close(3)
    
13.         call arburg(x,a,ef,eb,n,ip,ep,ierror)
    
14.         if(ierror.eq.0) goto 30
    
15.         write(*,*)'   stop at rotine aryuwa,    ierror=',ierror
    
16.         stop
    
17. 30      write(*,*)'   white noise variance=',ep
    
18.         write(*,*)'              k               a(k)'
    
19.         do 40 k=0,ip-1
    
20. 40      write(*,41)k,a(k)
    
21. 41      format('      ',i10,2f14.6)
    
22.         call ar1psd(a,ip,mfre,ep,ts)
    
23.         stop
    
24.         end

复制代码

子程序：

1.         subroutine arburg(x,a,ef,eb,n,ip,ep,ierror)
   
2. c----------------------------------------------------------------------
   
3. c   Routine ARBURG: To estimate the AR parameters by Burg algorithm.
   
4. c   Input Parameters:
   
5. c          N  : Number of data samples;
   
6. c          IP : Order of autoregressive model;
   
7. c          X  : Array of complex data samples, X(0) through X(N-1);
   
8. c   Output Parameters:
   
9. c         EP  : Real variable representing driving noise variance;
   
10. c          A  : Array of complex AR parameters A(0) to A(IP);
    
11. c    IERROR=0 : No error
    
12. c          =1 : Ep<=0 .
    
13. c                                      in chapter 12
    
14. c----------------------------------------------------------------------
    
15.         complex x(0:n-1),a(0:ip),ef(0:n-1),eb(0:n-1),sum
    
16.         ierror=1
    
17.         a(0)=1.
    
18.         r0=0.
    
19.         do 10 i=0,n-1
    
20.            r0=r0+x(i)*conjg(x(i))
    
21.            ef(i)=x(i)
    
22.            eb(i)=x(i)
    
23. 10      continue
    
24.         den=0.0
    
25.         r0=r0/float(n)
    
26.         sum=0.0
    
27.         do 20 i=1,n-1
    
28.            den=den+ef(i)*conjg(ef(i))+eb(i-1)*conjg(eb(i-1))
    
29.            sum=sum+ef(i)*conjg(eb(i-1))
    
30. 20      continue
    
31.         sum=-2.*sum
    
32.         a(1)=sum/den
    
33.         ep=r0*(1.-a(1)*conjg(a(1)))
    
34.         do 30 i=n-1,1,-1
    
35.            sum=ef(i)
    
36.            ef(i)=sum+a(1)*eb(i-1)
    
37.            eb(i)=eb(i-1)+conjg(a(1))*sum
    
38. 30      continue
    
39. c--------------------------------------------------------------------
    
40.         do 80 m=2,ip
    
41.            sum=0.0
    
42.            do 40 i=m,n-1
    
43.               sum=sum+ef(i)*conjg(eb(i-1))
    
44. 40         continue
    
45.            sum=-2.*sum
    
46.            den=(1.-a(m-1)*conjg(a(m-1)))*den-ef(m-1)*conjg(ef(m-1))
    
47.      *         -eb(n-1)*conjg(eb(n-1))
    
48.            a(m)=sum/den
    
49.            ep=ep*(1.-a(m)*conjg(a(m)))
    
50.            if(ep.le.0) return
    
51.            do 50 i=1,m-1
    
52.               x(i)=a(i)+a(m)*conjg(a(m-i))
    
53. 50         continue
    
54.            do 60 i=1,m-1
    
55.               a(i)=x(i)
    
56. 60         continue
    
57.            do 70 i=n-1,m,-1
    
58.               sum=ef(i)
    
59.               ef(i)=sum+a(m)*eb(i-1)
    
60.               eb(i)=eb(i-1)+conjg(a(m))*sum
    
61. 70         continue
    
62. 80      continue
    
63.         ierror=0
    
64.         return
    
65.         end

复制代码

halibut

第一次算变量den的值用公式
den=den+ef(i)*conjg(ef(i))+eb(i-1)*conjg(eb(i-1))
后面的用式子：
den=(1.-a(m-1)*conjg(a(m-1)))*den-ef(m-1)*conjg(ef(m-1))
     *         -eb(n-1)*conjg(eb(n-1))
这个式子是怎么推出来的？