由AR模型参数得到功率谱

VibInfo

1.       subroutine ar1psd(a,ip,mfre,ep,ts)
   
2. c-----------------------------------------------------------------------
   
3. c   Routine AR1PSD: To compute the power spectum by AR-model parameters.
   
4. c   Input parameters:
   
5. c          IP   : AR model order (integer)
   
6. c          EP   : White noise variance of model input (real)
   
7. c          Ts   : Sampling interval in seconds (real)
   
8. c          A    : Complex array of AR  parameters, A(0) to A(IP)
   
9. c   Output parameters:
   
10. c          PSDR : Real array of power spectral density values
    
11. c          PSDI : Real work array
    
12. c                                       in chapter 12
    
13. c-----------------------------------------------------------------------
    
14.         complex a(0:ip)
    
15.         dimension psdr(0:4095),psdi(0:4095)
    
16.         do 10 k=0,ip
    
17.            psdr(k)=real(a(k))
    
18.            psdi(k)=aimag(a(k))
    
19. 10      continue
    
20.         do 20 k=ip+1,mfre-1
    
21.            psdr(k)=0.
    
22.            psdi(k)=0.
    
23. 20      continue
    
24.         call relfft(psdr,psdi,mfre,-1)
    
25.         do 30 k=0,mfre-1
    
26.            p=psdr(k)**2+psdi(k)**2
    
27.            psdr(k)=ep*ts/p
    
28. 30      continue
    
29.         call psplot(psdr,psdi,mfre,ts)
    
30.         return
    
31.         end

复制代码

dxjuan115

PS:

relfft(psdr,psdi,mfre,-1)
psplot(psdr,psdi,mfre,ts)
是两个子程序吧，请问我可以在哪里找到？

风花雪月

1.         subroutine relfft(xr,xi,n,isign)
   
2. c----------------------------------------------------------------------
   
3. c  routinue RELFFT: To perform  the split-radix DIF fft algorithm.
   
4. c  input parameters:
   
5. c   xr,xi:real and image part of complex data for DFT/IDFT,n=0,...,N-1
   
6. c      n : the dimension of xr and xi ;
   
7. cISIGN=-1: For Forward Transform;
   
8. c     =+1: For Inverse Transform.
   
9. c  output parameters:
   
10. c   xr,xi:real and image part  of DFT/IDFT's result,n=0,...,N-1
    
11. c  Note  : N  must be a power of 2 .
    
12. c  From Ref. [7] of chapter 5, Modified by Lao Changan
    
13. c                                       in chapter 5
    
14. c----------------------------------------------------------------------
    
15.         dimension xr(0:n-1),xi(0:n-1)
    
16.         do 5 m=1,16
    
17.            n2=2**m
    
18.            if(n.eq.n2) go to 8
    
19. 5       continue
    
20.         write(*,*)'  N is not a power of 2'
    
21.         return
    
22. 8       n2=n*2
    
23.         es=-isign*atan(1.0)*8.0
    
24.         do 10 k=1,m-1
    
25.            n2=n2/2
    
26.            n4=n2/4
    
27.            e=es/n2
    
28.            a=0
    
29.            do 20 j=0,n4-1
    
30.               a3=3*a
    
31.               cc1=cos(a)
    
32.               ss1=sin(a)
    
33.               cc3=cos(a3)
    
34.               ss3=sin(a3)
    
35.               a=(j+1)*e
    
36.               is=j
    
37.               id=2*n2
    
38. 40            do 30 i0=is,n-1,id
    
39.               i1=i0+n4
    
40.               i2=i1+n4
    
41.               i3=i2+n4
    
42.               r1=xr(i0)-xr(i2)
    
43.               s1=xi(i0)-xi(i2)
    
44.               r2=xr(i1)-xr(i3)
    
45.               s2=xi(i1)-xi(i3)
    
46.               xr(i0)=xr(i0)+xr(i2)
    
47.               xi(i0)=xi(i0)+xi(i2)
    
48.               xr(i1)=xr(i1)+xr(i3)
    
49.               xi(i1)=xi(i1)+xi(i3)
    
50.               if(isign.eq.1)goto 41
    
51.               s3=r1-s2
    
52.               r1=r1+s2
    
53.               s2=r2-s1
    
54.               r2=r2+s1
    
55.               goto 42
    
56. 41            s3=r1+s2
    
57.               r1=r1-s2
    
58.               s2=-r2-s1
    
59.               r2=-r2+s1
    
60. 42            xr(i2)=r1*cc1-s2*ss1
    
61.               xi(i2)=-s2*cc1-r1*ss1
    
62.               xr(i3)=s3*cc3+r2*ss3
    
63.               xi(i3)=r2*cc3-s3*ss3
    
64. 30            continue
    
65.               is=2*id-n2+j
    
66.               id=4*id
    
67.               if(is.lt.n-1) goto 40
    
68. 20         continue
    
69. 10      continue
    
70. c
    
71. c       ------------ special last stage -------------------------
    
72.         is=0
    
73.         id=4
    
74. 50      do 60 i0=is,n-1,id
    
75.            i1=i0+1
    
76.            xtr=xr(i0)
    
77.            xti=xi(i0)
    
78.            xr(i0)=xtr+xr(i1)
    
79.            xi(i0)=xti+xi(i1)
    
80.            xr(i1)=xtr-xr(i1)
    
81.            xi(i1)=xti-xi(i1)
    
82. 60      continue
    
83.            is=2*id-2
    
84.            id=4*id
    
85.            if(is.lt.n-1) goto 50
    
86. 
    
87. 100     j=1
    
88.         n1=n-1
    
89.         do 104 i=1,n1
    
90.            if(i.ge.j) goto 101
    
91.            xtr=xr(j-1)
    
92.            xti=xi(j-1)
    
93.            xr(j-1)=xr(i-1)
    
94.            xi(j-1)=xi(i-1)
    
95.            xr(i-1)=xtr
    
96.            xi(i-1)=xti
    
97. 101        k=n/2
    
98. 102        if(k.ge.j) goto 103
    
99.                 j=j-k
    
100.                 k=k/2
     
101.                 goto 102
     
102. 103        j=j+k
     
103. 104     continue
     
104.         if(isign.eq.-1)return
     
105.         do 105 i=0,n-1
     
106.         xr(i)=xr(i)/n
     
107.         xi(i)=xi(i)/n
     
108. 105     continue
     
109.         return
     
110.         end

复制代码

风花雪月

1.       subroutine psplot(psdr,psdi,mfre,ts)
   
2. c-----------------------------------------------------------------------
   
3. c   Routine PSPLOT: To plot the normalized power spectum curve on the
   
4. c                      normalized frequency axis from -.5 to  +.5 .
   
5. c        mfre : Points in frequency axis and must be the power of 2.
   
6. c        ts    : Sampling interval in seconds (real).
   
7. c        psdr : Real array of power spectral density values.
   
8. c        psdi : Real work array.
   
9. c                                            in Chapter 11
   
10. c-----------------------------------------------------------------------
    
11.         dimension psdr(0:mfre-1),psdi(0:mfre-1)
    
12.         m2=mfre/2
    
13.         do 40 k=0,m2-1
    
14.            psdi(k)=psdr(k)
    
15.            psdr(k)=psdr(k+m2)
    
16.            psdr(k+m2)=psdi(k)
    
17. 40      continue
    
18.         pmax=psdr(0)
    
19.         do 50 k=1,mfre-1
    
20.            if(psdr(k).gt.pmax) pmax=psdr(k)
    
21. 50      continue
    
22.         do 60 k=0,mfre-1
    
23.            psdr(k)=psdr(k)/pmax
    
24.            if(psdr(k).le.0.0)psdr(k)=.000001
    
25. 60      continue
    
26.         fs=1./ts
    
27.         fs=fs/float(mfre)
    
28.         open(3,file='psd.dat',status='new')
    
29.         do 100 k=0,mfre-1
    
30.            faxis=fs*(k-m2)
    
31.            write(3,*)faxis,10.*alog10(psdr(k))
    
32. 100      continue
    
33.         close(3)
    
34.         return
    
35.         end

复制代码

风花雪月

原帖由 dxjuan115 于 2006-11-13 16:00 发表
PS:

relfft(psdr,psdi,mfre,-1)
psplot(psdr,psdi,mfre,ts)
是两个子程序吧，请问我可以在哪里找到？

你看一下是否是上面给的两个程序

lizhimin

收藏了，谢谢

zmxzhang

好东西，收藏，谢谢