第七章譜分析

1、第七章譜分析7快速傅立葉變換7.2.1功能用FFT計(jì)算離散傅立葉（Fourier）變換7.2.3子程序語句SUBROUTINE FFT(N,K,PR,PI, FR,FI,L,IL)7.2.4啞元說明N整型變量，輸入?yún)?shù)（即INTENT(IN)），序列的長度。K整型變量，輸入?yún)?shù)（即INTENT(IN)），滿足N=2。PR實(shí)型一維數(shù)組，長度為N，輸入兼輸出參數(shù)（即INTENT(INOUT)）。調(diào)用時(shí)存放N個(gè)輸入樣本的實(shí)部（L=0時(shí)）或傅立葉變換的N個(gè)實(shí)部（L=1時(shí)）；返回時(shí)存放傅立葉變換的模（L=0時(shí)）或逆傅立葉變換的模（L=1時(shí)）。PI實(shí)型一維數(shù)組，長度為N，輸入兼輸出參數(shù)（即INTENT(I

2、NOUT)）。調(diào)用時(shí)存放N個(gè)輸入樣本的虛部（L=0時(shí)）或傅立葉變換的N個(gè)虛部（L=1時(shí)）；返回時(shí)存放傅立葉變換的幅角（L=0時(shí)）或逆傅立葉變換的幅角（L=1時(shí)）。FR實(shí)型一維數(shù)組，長度為N，輸出參數(shù)。存放傅立葉變換（L=0時(shí)）或逆傅立葉變換（L=1時(shí)）的實(shí)部。FI實(shí)型一維數(shù)組，長度為N，輸出參數(shù)。存放傅立葉變換（L=0時(shí)）或逆傅立葉變換（L=1時(shí)）的虛部。L整型變量，輸入?yún)?shù)。如果L=0，則為計(jì)算傅立葉變換；如果L=1，則為計(jì)算逆傅立葉變換。IL整型變量，輸入?yún)?shù)。如果IL=0，為不計(jì)算傅立葉變換或逆變換的模與幅角；如果IL0為計(jì)算傅立葉變換或逆傅立葉變換的模與幅角。一、子程序（子程序名：FF

3、T）SUBROUTINE FFT(N,K,PR,PI,FR,FI,L,IL)REAL*8,DIMENSION(N): PR,PI,FR,FIREAL*8: P,Q,S,VR,VI,PODDR,PODDIREAL*8,PARAMETER:PAI=3.1415926DO IT=0,N-1M=ITIS=0DO I=0,K-1J=M/2IS=2*IS+(M-2*J)M=JEND DOFR(IT+1)=PR(IS+1)FI(IT+1)=PI(IS+1)END DOPR(1)=1.1 kPI(1)=0.PR(2)=COS(2*PAI/N)PI(2)=-SIN(2*PAI/N)IF(L.NE.0)PI(2)

4、=-PI(2)DO I=3,NP=PR(I-1)*PR(2)Q=PI(I-1)*PI(2)S=(PR(I-1)+PI(I-1)*(PR(2)+PI(2)PR(I)=P-QPI(I)=S-P-QEND DODO IT=0,N-2,2VR=FR(IT+1)VI=FI(IT+1)FR(IT+1)=VR+FR(IT+2)FI(IT+1)=VI+FI(IT+2)FR(IT+2)=VR-FR(IT+2)FI(IT+2)=VI-FI(IT+2)END DOM=N/2NV=2DO L0=K-2,0,-1M=M/2NV=2*NVDO IT=0,(M-1)*NV,NVDO J=0,(NV/2)-1P=PR(M*J

5、+1)*FR(IT+J+1+NV/2) Q=PI(M*J+1)*FI(IT+J+1+NV/2) S=PR(M*J+1)+PI(M*J+1)S=S*(FR(IT+J+1+NV/2)+FI(IT+J+1+NV/2) PODDR=P-QPODDI=S-P-QFR(IT+J+1+NV/2)=FR(IT+J+1)-PODDR FI(IT+J+1+NV/2)=FI(IT+J+1)-PODDI FR(IT+J+1)=FR(IT+J+1)+PODDR FI(IT+J+1)=FI(IT+J+1)+PODDIEND DOEND DOEND DOIF(L.NE.0)THENDO I=1,NFR(I)=FR(I)/N

6、FI(I)=FI(I)/NEND DOENDIFIF(IL.NE.0)THENDO I=1,NPR(I)=SQRT(FR(I)*FR(I)+FI(I)*FI(I)PI(I)=ATAN(FI(I)/FR(I)*180/PAIEND DOEND IFEND7.2.6例對(duì)某海區(qū)連續(xù)64個(gè)月的海表溫度作傅立葉變換。計(jì)算：（1）Pi的離散傅立葉變換Fi，且求出Fi的模與幅角。即L=0，IL=1。（2）Fi的逆傅立葉變換Pi的模。即L=1，IL=1。主程序?yàn)椋篜ROGRAM MAININTEGER,PARAMETER:K=6INTEGER,PARAMETER:N=2*KREAL*8,DIMENSION(N

7、): PR,PI,FR,FIOPEN(10,FILE='AA2.DAT')READ(10,'(F8.2)')PRCLOSE(10)PI=0.OPEN(12,FILE='FFT.DAT')WRITE(12,'(" 序列P(I)即原始資料")')WRITE(12,20)(PR(I),I=1,N)CALL FFT(N,K,PR,PI,FR,FI,0,1)WRITE(12,'(" 序列P(I)FFT變換F(I)的實(shí)部")')WRITE(12,20)(FR(I),I=1,N)WRITE

8、(12,'(" 序列P(I)FFT變換F(I)的虛部")')WRITE(12,20)(FI(I),I=1,N)WRITE(12,'(" 序列P(I)FFT變換F(I)的模")')WRITE(12,20)(PR(I),I=1,N)WRITE(12,'(" 序列P(I)FFT變換F(I)的幅角")')WRITE(12,20)(PI(I),I=1,N)CALL FFT(N,K,FR,FI,PR,PI,1,1)WRITE(12,'(" F(I)的FFT逆變換的模")&

9、#39;)WRITE(12,20)(FR(I),I=1,N)20 FORMAT(1X,8F8.2)CLOSE(12)END計(jì)算結(jié)果為：序列P(I)即原始資料:16.80 15.35 17.00 22.50 23.50 27.00 27.60 27.15 24.00 20.85 18.25 16.20 14.30 16.55 3 28.00 21.1024.00 26.25 27.80 27.30 27.05 25.50 23.80 19.95 15.60 17.00 19.70 20.90 24.00 24.80 26.95 26.70 27.40 24.85 22.20 18.90 15.8

10、0 13.55 17.60 21.75 25.00 26.20 26.95 27.00 26.35 24.60 21.55 17.85 15.60 18.05 18.90 21.90 24.35 26.20 26.80 26.90 28.05 25.60 22.00 17.80 16.20 15.20 17.60 20.00 序列P(I)FFT變換F(I)的實(shí)部:1412.15 -20.84 -23.73 -20.72 -44.19 -126.60 49.35 7.35 15.94 6.15 5.83 .82 -1.33 -1.54 -5.92 -.52 -.80 6.91 5.43-3.64

11、 4.88 .731.65 -6.68 1.75-3.64 .55 3.36-.80 -.52 -5.9215.94 7.35 49.35序列P(I)FFT變換F(I)的虛部:.00 -2.68 3.874.15 12.41 -.998.35 4.08 2.53-1.95 1.69 -3.61.00 -.34 2.301.95 -.25 -3.48-8.35 -1.72 -15.12-4.15 -22.12 -53.49序列P(I)FFT變換F(I)的模:1412.15 21.02 24.0416.47 13.85 5.928.39 8.02 5.994.13 5.17 3.681.65 6.

12、69 2.894.13 .60 4.848.39 1.80 16.2416.47 23.31 72.78序列P(I)FFT變換F(I)的幅角:.00 7.32 -9.2614.58 63.64 -9.67-84.53 30.58 24.9928.20 19.10 -78.61.00 2.95 52.73-28.20 -24.05 -46.0584.53 73.08 68.62-14.58 -71.62 -47.30F(I)的FFT逆變換的模:16.80 15.35 17.00 -9.17 1.02 -1.72 3.36 .55 -3.54 -3.50 -4.93 1.75 -6.68 -4.9

13、3 -3.50 -3.54 .73 4.88 -1.72 1.02 -9.17 5.43 6.91 -1.54 -1.33 .82 5.83 6.15 -126.60 -44.19 -20.72 -23.73 -20.84 -18.88 -34.19 -97.58 53.49 22.12 37.24 6.40 1.41 15.12 1.72 8.11 4.01 3.62 3.48 .25 .22 -4.79 -6.23 -2.30 .34 6.23 4.79 -.22 3.61 -1.69 -3.62 -4.01 -8.11 -2.53 -4.08 -1.41 -6.40 -37.24 .99

14、 -12.41 97.58 34.19 18.88 -3.87 2.68 28.04 55.88 159.84 72.78 23.31 37.25 6.54 2.09 16.24 1.80 12.24 4.14 4.01 4.84 .60 3.54 5.93 7.94 2.89 6.69 7.94 5.93 3.54 3.68 5.17 4.01 4.14 12.24 5.99 8.02 2.09 6.54 37.25 5.92 13.85 159.84 55.88 28.04 24.04 21.02 42.34 37.73 37.62 47.30 71.62 88.74 -78.24 -42

15、.43 -68.62 -73.08 -41.51 75.71 -64.61 46.05 24.05 -3.60 53.87 51.65 -52.73 -2.95 -51.65 -53.87 3.60 78.61 -19.10 64.61 -75.71 41.51 -24.99 -30.58 42.43 78.24 -88.74 9.67 -63.64 -37.62 -37.73 -42.34 9.26 -7.32 22.50 23.50 27.00 27.60 28.0027.15 24.00 20.85 18.25 16.20 14.30 16.55 21.10 24.00 26.25 27

16、.80 27.30 27.05 25.50 23.80 19.95 15.60 17.00 19.70 20.90 24.00 24.80 26.95 26.70 27.40 24.85 22.20 18.90 15.80 13.55 17.60 21.75 25.00 26.20 26.95 27.00 26.35 24.60 21.55 17.85 15.60 18.05 18.90 21.90 24.35 26.20 26.80 26.90 28.05 25.60 22.00 17.80 16.20 15.20 17.60 20.00 由計(jì)算結(jié)果可見，逆變換所得到的值和原輸入的資料完全相

17、同。功率譜分析7.3.1功能計(jì)算序列的平滑功率譜7.3.2方法說明7.3.3子程序語句SUBROUTINE SPECTRUM1(X,N,M,R,S)7.3.4啞元說明N整型變量，輸入?yún)?shù)。序列的長度。M整型變量，輸入?yún)?shù)。最大落后時(shí)間長度。X實(shí)型一維數(shù)組，長度為N，輸入?yún)?shù)，存放觀測的序列值。R實(shí)型一維數(shù)組，長度為M+1，輸出參數(shù)，存放M+1個(gè)自相關(guān)系數(shù)。 S實(shí)型一維數(shù)組，長度為M+1，輸出參數(shù)，存放M+1個(gè)功率譜值。7.3.5子程序（子程序名為：SPECTRUM1）SUBROUTINE SPECTRUM1(X,N,M,R,S) REAL(4),DIMENSION(N):X REAL(4),D

18、IMENSION(0:M):R,B,S REAL(4): XBAR INTEGER: TAO,T REAL(4),PARAMETER:PI=3.1415926 XBAR=0 DO I=1,N XBAR=XBAR+X(I) END DO XBAR=XBAR/N S1=0 DO I=1,N S1=S1+(X(I)-XBAR)*2 END DO DO TAO=0,M R(TAO)=0 DO T=1,N-TAO R(TAO)=R(TAO)+(X(T)-XBAR)*(X(T+TAO)-XBAR)/S1 END DO R(TAO)=R(TAO)/(N-TAO) 5END DO DO L=0,M IF(L.

19、EQ.0.OR.L.EQ.M)THEN B(L)=1 ELSE B(L)=1./2 END IF S(L)=R(0) DO TAO=1,M-1 S(L)=S(L)+R(TAO)*(1+COS(PI*TAO/M)*COS(L*PI*TAO/M) END DO S(L)=B(L)*S(L)/M END DO END7.3.6例計(jì)算某海區(qū)1971年到1982年逐月的海表溫度的月平均值的功率譜。序列長度為N=12*12=144，取M=18。海表溫度值存放在文件SST.DAT中。PROGRAM MAININTEGER,PARAMETER(N=144,M=18)REAL(4),DIMENSION(N):X

20、REAL(4),DIMENSION(0:M):R,SOPEN(10,FILE='SST.DAT')READ(10,'(12F8.2)')XCLOSE(10)CALL SPECTRUM1(X,N,M,R,S)OPEN(12,FILE=' SPECTRUM1.DAT')DO I=0,MWRITE(12, (2F8.3)R(I),S(I)END DOEND附： 海表溫度(SST.DAT)序列為：16.80 15.35 17.00 22.50 23.50 27.00 27.60 28.00 27.15 24.00 20.85 18.25 16.20 14

21、.30 16.55 21.10 24.00 26.25 27.80 27.30 27.05 25.50 23.80 19.95 15.60 17.00 19.70 20.90 24.00 24.80 26.95 26.70 27.40 24.85 22.20 18.90 15.80 13.55 17.60 21.75 25.00 26.20 26.95 27.00 26.35 24.60 21.55 17.85 15.60 18.05 18.90 21.90 24.35 26.20 26.80 26.90 28.05 25.60 22.00 17.80 16.20 15.20 17.60 20

22、.00 23.75 25.20 27.00 27.80 26.90 24.40 21.00 17.80 14.00 13.55 19.95 23.00 25.15 26.80 27.00 27.10 26.80 25.50 22.20 19.50 18.00 17.80 18.95 21.70 23.40 27.35 28.00 27.80 27.20 25.00 22.20 19.95 18.95 19.00 20.50 22.20 23.35 25.55 27.90 27.80 28.00 24.60 22.50 19.20 17.70 15.10 16.50 22.00 24.00 28

23、.00 28.60 27.90 27.00 25.40 23.00 21.30 18.50 18.00 19.00 23.25 24.25 25.40 28.10 28.50 26.70 25.70 22.00 18.00 18.00 17.00 18.00 20.00 24.05 25.50 27.55 27.50 26.60 26.00 23.50 20.00 輸出結(jié)果為：L R(L) S(L)0 1.0000 .00501 .8353 .00902 .4637 .23823 -.0105 .46644 -.4620 .23505 -.7870 .00876 -.9002 .01207 -

24、.7887 .00798 -.4704 .00329 -.0155 .002410 .4533 .002111 .8014 .002312 .9361 .001913 .8085 .001014 .4561 .000915 .0019 .001016 -.4496 .001117 -.7664 .001418 -.8848 .000874 交叉譜分析7.4.1功能計(jì)算兩個(gè)不同的時(shí)間序列x1(t)，x2 (t)的交叉譜。7.4.2方法說明為不同的復(fù)譜，故交叉譜仍是復(fù)譜，它可以寫成實(shí)部與虛部之和的形式，即：7.4.3子程序語句SUBROUTINE SPECTRUM2(N,M,X1,X2, R12,

25、R21,P11,P22,P12,Q12,R212,THITA)7.4.4啞元說明N整型變量，輸入?yún)?shù)，序列長度。M整型變量，輸入?yún)?shù)，最大落后長度。X1實(shí)型一維數(shù)組，輸入?yún)?shù)，長度為N，存放第一要素觀測序列值。X2實(shí)型一維數(shù)組，輸入?yún)?shù)，長度為N，存放第二要素觀測序列值。R11實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X1的自相關(guān)系數(shù)。R22實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X2的自相關(guān)系數(shù)。R12實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放落后相關(guān)系數(shù)（X2落后于X1）。 R21實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放落后相關(guān)系數(shù)（X1落后于X2）。 P11實(shí)型一維數(shù)組，輸出參數(shù)，長

26、度為M+1，存放X1的功率譜。P22實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X2的功率譜。P12實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X1和X2的協(xié)譜。Q12實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X1和X2的正交譜。 R212實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X1和X2的凝聚譜。 THITA實(shí)型一維數(shù)組，輸出參數(shù)，長度為M+1，存放X1和X2的位相差譜。7.4.5子程序（子程序名為：SPECTRUM2）SUBROUTINE SPECTRUM2(N,M,X1,X2, R12,R21,P11,P22,P12,Q12,R212,THITA) REAL(4),DIMENSION(N)

27、:X1,X2 REAL(4),DIMENSION(0:M):R12,R21,R11,R22,P12,Q12,P22,P11,R212,THITA,B REAL(4): X1BAR,X2BAR INTEGER: TAO,T REAL(4),PARAMETER:PI=3.1415926 X1BAR=0 X2BAR=0 DO I=1,N X1BAR=X1BAR+X1(I) X2BAR=X2BAR+X2(I) END DO X1BAR=X1BAR/N X2BAR=X2BAR/N S1=0 S2=0 DO I=1,N S1=S1+(X1(I)-X1BAR)*2 S2=S2+(X2(I)-X2BAR)*2

28、 END DO S1=SQRT(S1/N) S2=SQRT(S2/N) DO TAO=0,M R11(TAO)=0 R22(TAO)=0 R12(TAO)=0 R21(TAO)=0 DO T=1,N-TAO R11(TAO)=R11(TAO)+(X1(T)-X1BAR)/S1*(X1(T+TAO)-X1BAR)/S1 R22(TAO)=R22(TAO)+(X2(T)-X2BAR)/S2*(X2(T+TAO)-X2BAR)/S2 R12(TAO)=R12(TAO)+(X1(T)-X1BAR)/S1*(X2(T+TAO)-X2BAR)/S2 R21(TAO)=R21(TAO)+(X2(T)-X2B

29、AR)/S2*(X1(T+TAO)-X1BAR)/S1 END DO R12(TAO)=R12(TAO)/(N-TAO) R21(TAO)=R21(TAO)/(N-TAO) END DO DO L=0,M IF(L.EQ.0.OR.L.EQ.M)THEN B(L)=0.5 ELSE B(L)=1. END IF P11(L)=R11(0) P22(L)=R22(0) 8P12(L)=R12(0) Q12(L)=0 DO TAO=1,M-1 P11(L)=P11(L)+R11(TAO)*(1+COS(PI*TAO/M)*COS(L*PI*TAO/M) P22(L)=P22(L)+R22(TAO)

30、*(1+COS(PI*TAO/M)*COS(L*PI*TAO/M) P12(L)=P12(L)+0.5*(R12(TAO)+R21(TAO)*(1+COS(PI*TAO/M)*COS(L*PI*TAO/M) Q12(L)=Q12(L)+0.5*(1+COS(PI*TAO/M)*SIN(PI*L*TAO/M)*(R12(TAO)-R21(TAO) END DO P11(L)=B(L)*P11(L)/M P22(L)=B(L)*P22(L)/M P12(L)=B(L)*P12(L)/M Q12(L)=B(L)*Q12(L)/M END DO DO L=1,M-1 R212(L)=(P12(L)*2

31、+Q12(L)*2)/(P11(L)*P22(L) THITA(L)=180/PI*ATAN(Q12(L)/P12(L) END DO END7.4.6例計(jì)算某海區(qū)10年逐月的海表溫度和其水面上的氣溫的月平均值的交叉譜。序列長度為N=12*10=120，取M=12。海表溫度值存放在文件aa2.dat中。PROGRAM MAINPARAMETER(N=120,M=12)REAL(4),DIMENSION(N):X1,X2REAL(4),DIMENSION(0:M):R12,R21,P11,P22,P12,Q12,R212,THITAOPEN(10,FILE='AA2.DAT')D

32、O I=1,NREAD(10,'(2F8.2)')X1(I),X2(I)END DOCLOSE(10)CALL SPECTRUM2(X1,X2,N,M,R12,R21,P11,P22,P12,Q12,R212,THITA)OPEN(12,FILE=' SPECTRUM2.DAT')WRITE(12,'(" P11 P22 P12 Q12 R212 THITA ")') DO I=0,MWRITE(12, (6F10.5)P11(I),P22(I),P12(I),Q12(I),R212(I),THITA(I)END DOEND計(jì)

33、算結(jié)果:P11 P22 P12 Q12 R212 THITA1.51247 .96051 .00389 .00000 .00000 .0000028.94028 28.65097 .23889 -.01301 .00007 -3.1184154.76953 55.47872 .47035 -.02624 .00007 -3.1926528.95655 29.08382 .24203 -.01254 .00007 -2.966162.52968 1.94579 .01212 .00164 .00003 7.710211.25791 .96728 .00703 .00145 .00004 11.6

34、7822.55827 .56115 .00248 .00080 .00002 17.81564.40535 .52044 .00187 .00020 .00002 6.15760.28200 .49032 .00138 -.00036 .00001 -14.69649.18011 .38023 .00066 -.00069 .00001 -46.52385.21489 .37978 .00105 -.00050 .00002 -25.59904.25977 .39727 .00135 .00017 .00002 6.97895.13311 .18373 .00055 .00000 .00000

35、 .000007最大熵譜分析7.5.1功能用最大熵譜法作譜估計(jì)。7.5.2方法說明7.5.3子程序語句SHANGPU(X,N,M,K0,SE,T,OMEGA,SIGMA2,FAI,B)7.5.4啞元說明X實(shí)型一維數(shù)組，輸入?yún)?shù)，長度為N，存放要素觀測序列值。N整型變量，輸入?yún)?shù)，序列長度。M整型變量，輸入?yún)?shù)，最大落后長度。K0輸出參數(shù)，存放最佳的階數(shù)。SE實(shí)型一維數(shù)組，輸出參數(shù)，長度為M，存放熵譜值。T實(shí)型一維數(shù)組，輸出參數(shù)，長度為M，存放熵譜值對(duì)應(yīng)的周期。SIGMA2實(shí)型一維數(shù)組，輸出參數(shù)，長度為0:N-1，存放誤差的方差。FAI實(shí)型一維數(shù)組，輸出參數(shù)，長度為N-1，存放Akaike的“最

36、終預(yù)報(bào)誤差”。 B實(shí)型二維數(shù)組，輸出參數(shù)，大小為（N-1，N-1），存放遞推的回歸系數(shù)。7.5.5子程序(子程序名為：SHANGPU)SUBROUTINE SHANGPU(X,N,M,K0,SE,T,OMEGA,SIGMA2,FAI,B) INTEGER:N,M REAL(8),DIMENSION(N):X REAL(8),DIMENSION(N-1,N-1):B REAL(8),DIMENSION(0:N-1):SIGMA2 REAL(8),DIMENSION(N-1):FAI REAL(8),DIMENSION(0:M):SE,T,OMEGA REAL(8):A1,A2,A3,A4 REA

37、L(8):PI PI=4*ATAN(1.) B=0 A3=0 DO I=2,N A3=A3+X(I)*X(I-1) END DO A4=(X(1)*X(1)+X(N)*X(N)/2 DO I=2,N-110A4=A4+X(I)*2 END DO B(1,1)=A3/A4 DO K=1,N-2 A1=0 A2=0 DO I=K+2,N A3=X(I) A4=X(I-K-1) DO J=1,K A3=A3-B(K,J)*X(I-J) A4=A4-B(K,J)*X(I-K-1+J) END DO A1=A1+A3*A4 A2=A2+A3*2+A4*2 END DO B(K+1,K+1)=2*A1/A2 DO J=1,K B(K+1,J)=B(K,J)-B(K+1,K+1)*B(K,K+1-J) END DO END DO SIGMA2(0)=0 DO I=1,N SIGMA2(0)=SIGMA2(0)+X(I)*X(I) END DO SIGMA2(0)=SIGMA2(0)/N DO K=1,N-1 SIGMA2(K)=SIGMA2(K-1)*(1-B(K,K)*2) FAI(K)=SIGMA2(K)*(N+K)/(N-K) END D