7．4 交叉谱分析 7.4.1功能 计算两个不同的时间序列x1(t)，x2 (t)的交叉谱。 7.4.2方法说明 为不同的复谱，故交叉谱仍是复谱，它可以写成实部与虚部之和的形式，即： 7.4.3子程序语句 SUBROUTINE SPECTRUM2(N,M,X1,X2, R12,R21,P11,P22,P12,Q12,R212,THITA) 7.4.4哑元说明 N——整型变量，输入参数，序列长度。 M——整型变量，输入参数，最大落后长度。 X1——实型一维数组，输入参数，长度为N，存放第一要素观测序列值。 X2——实型一维数组，输入参数，长度为N，存放第二要素观测序列值。 R11——实型一维数组，输出参数，长度为M+1，存放X1的自相关系数。 R22——实型一维数组，输出参数，长度为M+1，存放X2的自相关系数。 R12——实型一维数组，输出参数，长度为M+1，存放落后相关系数（X2落后于X1）。 R21——实型一维数组，输出参数，长度为M+1，存放落后相关系数（X1落后于X2）。 P11——实型一维数组，输出参数，长度为M+1，存放X1的功率谱。 P22——实型一维数组，输出参数，长度为M+1，存放X2的功率谱。 P12——实型一维数组，输出参数，长度为M+1，存放X1和X2的协谱。 Q12——实型一维数组，输出参数，长度为M+1，存放X1和X2的正交谱。 R212——实型一维数组，输出参数，长度为M+1，存放X1和X2的凝聚谱。 THITA——实型一维数组，输出参数，长度为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)::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 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)-X2BAR)/S2*(X1(T+TAO)-X1BAR)/S1 END DO R11(TAO)=R11(TAO)/(N-TAO) !wangzhe R22(TAO)=R22(TAO)/(N-TAO) !wangzhe 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) P12(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)*(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+Q12(L)**2)/(P11(L)*P22(L)) THITA(L)=180/PI*ATAN(Q12(L)/P12(L)) END DO END 7.4.6例 计算某海区10年逐月的海表温度和其水面上的气温的月平均值的交叉谱。序列长度为N=12*10=120，取M=12。海表温度值存放在文件aa2.dat中。 PROGRAM MAIN PARAMETER(N=120,M=12) REAL(4),DIMENSION(N)::X1,X2 REAL(4),DIMENSION(0:M)::R12,R21,P11,P22,P12,Q12,R212,THITA OPEN(10,FILE='AA2.DAT') DO I=1,N READ(10,'(2F8.2)')X1(I),X2(I) END DO CLOSE(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,M WRITE(12,’ (6F10.5)’)P11(I),P22(I),P12(I),Q12(I),R212(I),THITA(I) END DO END 计算结果: P11 P22 P12 Q12 R212 THITA 1.51247 .96051 .00389 .00000 .00000 .00000 28.94028 28.65097 .23889 -.01301 .00007 -3.11841 54.76953 55.47872 .47035 -.02624 .00007 -3.19265 28.95655 29.08382 .24203 -.01254 .00007 -2.96616 2.52968 1.94579 .01212 .00164 .00003 7.71021 1.25791 .96728 .00703 .00145 .00004 11.67822 .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 .00000 本文来源：https://www.wddqw.com/doc/b011e2e433d4b14e852468ca.html