三角级数法模拟伪随机序列求助

行者无疆TJ

以下是自编的在0-1间伪随机数的基础上利用三角级数法模拟的随机序列.已知谱密度为Sx(w)=0.09
50.24<w<62.8即8<f<10.0
C-------------------------------------------------------------------------
C           PROGRAM: 白噪声随机数的产生
C-------------------------------------------------------------------------
        PARAMETER (NL1=1000000,NL2=2*NL1)
        REAL*8 XA(NL1),F(NL1)
        REAL*8 S(NL2),W(NL2),A(NL2)
        REAL*8 DELTAF,P,Y,U,V,FC
        WRITE(*,9997)
9997 FORMAT('IDUM=')
      READ(*,*) IDUM
        WRITE(*,9996)
9996 FORMAT('FC=')
      READ(*,*)  FC
      WRITE(*,9990)
9990 FORMAT('FL=')
      READ(*,*)  FL
        WRITE(*,9995)
9995 FORMAT('N=')
      READ(*,*)N
        U=8.0D0*ATAN(1.0D0)
      DELTAF=FC/N
        DO 10 I=1,N
        A(I)=RAN1(IDUM)
        WRITE(*,9994) A(I)
   10 CONTINUE
9994 FORMAT(F12.6)
        DO 20 R=1,2*N
      XA(R)=0.0D0
        DO 30 K=1,N-1
        F(K)=K*DELTAF
        IF (F(K).GE.FL) THEN
        IF (F(K).LE.FC)  THEN
        S(K)=(9.0D-2)*U
        ELSE
         S(K)=0.0D0
        END IF
        END IF
        Y=DSQRT(4.0D0*S(K)*DELTAF)
        P=DCOS(U*K*(R-1)/(2.0D0*N)+A(K)*U)
        XA(R)=XA(R)+Y*P
   30 CONTINUE
   20 CONTINUE
      OPEN(2,FILE='BZSHPL3.DAT',STATUS='NEW',
     *ACCESS='DIRECT',FORM='FORMATTED',RECL=15)
      DO 40  K=1,N
        WRITE(2,9992,REC=K) F(K)
   40 CONTINUE
      CLOSE(2)
9992 FORMAT(F10.6)
      OPEN(3,FILE='BZSSJS3.DAT',STATUS='NEW',
     *ACCESS='DIRECT',FORM='FORMATTED',RECL=15)
        DO 50  I=1,2*N
      WRITE(3,9991,REC=I)XA(I)
9991 FORMAT(F10.6)
   50 CONTINUE
      CLOSE(3)
      END
C
C       
C   
      FUNCTION RAN1(IDUM)
        INTEGER IDUM
        INTEGER M1,M2,M3,IA1,IA2,IA3,IC1,IC2,IC3,IFF,IX1,IX2,IX3
        REAL*8 R(97),RM1,RM2,RAN1
        PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.85802D-6)
        PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.43738D-6)
        PARAMETER (M3=243000,IA3=4561,IC3=51349)
        DATA IFF/0/
C
        IF (IDUM.LT.0 .OR. IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IA2*IX2+IC2,M2)
        R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1
11   CONTINUE
        IDUM=1
        ENDIF
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IA2*IX2+IC2,M2)
        IX3=MOD(IA3*IX3+IC3,M3)
        J=1+(97*IX3)/M3
        IF (J.GT.97 .OR. J.LT.1) PAUSE
        RAN1=R(J)
        R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1
        RETURN
        END

在此基础上进行付立叶计算,总不能返回到原谱密度?

风花雪月

在此基础上进行付立叶计算,总不能返回到原谱密度?

能够说的具体一点，什么地方要返回原谱密度？

行者无疆TJ

呵呵,总算有点回应,谢谢,表达能力实需提高了
C-------------------------------------------------------------------------
C           PROGRAM: 白噪声随机序列的产生
C-------------------------------------------------------------------------
       　　 PARAMETER (NL1=1000000,NL2=2*NL1)　　　　　　　　％定义数据类型
       　　  REAL*8 XA(NL1),F(NL1)
     　　   REAL*8 S(NL2),W(NL2),A(NL2)
       　　 REAL*8 DELTAF,P,Y,U,V,FC
      　　  WRITE(*,9997)                                           ％读入数据IDUM,N,FC,IDUM输入为一，主要
9997  　　FORMAT('IDUM=')                                           后面的函数获取０－１间同一随机数而备，
                READ(*,*) IDUM
        　　　WRITE(*,9996)
9996        FORMAT('FC=')
                READ(*,*)  FC                                               FC为最高频率
     　　   WRITE(*,9995)
9995 　　FORMAT('N=')
     　　　 READ(*,*)N　　　　　　　　　　　　　　　　Ｎ将频率分为Ｎ段
      　　  U=8.0D0*DATAN(1.0D0) 　　　　　　　　　　　Ｕ即为２pi
               DELTAF=FC/N
     　　   DO 10 I=1,N                                                      %从函数RAN1调取０－１间同一随机数
       　　 A(I)=RAN1(IDUM)
       　　 WRITE(*,9994) A(I)
   10         CONTINUE
9994        FORMAT(F12.6)
      　　　  DO 20 R=1,2*N                                             %20循环：获取一系列时间序列
               　 XA(R)=0.0D0
      　　　  DO 30 K=0,N-1                                              ％30循环：每一循环获取一时间序列
     　　　   F(K)=K*DELTAF
    　   　　 IF (F(K+1).GE.8.0D0) THEN
     　　　   IF (F(K+1).LE.1.0D1)  THEN
       　　　 S(K+1)=9.0D-2*U
      　　　  ELSE
     　　   　 S(K+1)=0.0D0
     　　 　  END IF
     　　 　  END IF
      　　  　Y=DSQRT(4.0D0*S(K+1)*DELTAF)
      　　 　 P=DCOS(U*K*(R-1)/(2.0D0*N)+A(K+1)*U)
       　　　 XA(R)=XA(R)+Y*P
   30         CONTINUE
   20          CONTINUE
                OPEN(3,FILE='BZSSJS3.DAT',STATUS='NEW',                   %随机序列以文件输出
       　     *ACCESS='DIRECT',FORM='FORMATTED',RECL=15)
      　　　  DO 50  I=1,2*N
             WRITE(3,9991,REC=I)XA(I)
9991      FORMAT(F10.6)
   50        CONTINUE
               CLOSE(3)
                END
C
C　　以下为获取０－１间同一随机数的函数        （现成的，应该没问题）
C   
    　　　  FUNCTION RAN1(IDUM)
        INTEGER IDUM
        INTEGER M1,M2,M3,IA1,IA2,IA3,IC1,IC2,IC3,IFF,IX1,IX2,IX3
        REAL*8 R(97),RM1,RM2,RAN1
        PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.85802D-6)
        PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.43738D-6)
        PARAMETER (M3=243000,IA3=4561,IC3=51349)
        DATA IFF/0/
C
        IF (IDUM.LT.0 .OR. IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IA2*IX2+IC2,M2)
        R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1
11   　CONTINUE
        IDUM=1
        ENDIF
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IA2*IX2+IC2,M2)
        IX3=MOD(IA3*IX3+IC3,M3)
        J=1+(97*IX3)/M3
        IF (J.GT.97 .OR. J.LT.1) PAUSE
        RAN1=R(J)
        R(J)=(DFLOAT(IX1)+DFLOAT(IX2)*RM2)*RM1
        RETURN
　        END
这个程序只是想通过星谷胜的三角级数模拟出随机序列,反算谱是另一个程序.想要模拟的谱如图１,谱密度Sx(w)=0.09  (50.24<w<62.8)
书中直接ＤＦＴ变换模拟结果为图２，　Sx(f)在1.1左右，如何理解？而我模拟出来的谱值１１左右？

[ 本帖最后由 行者无疆TJ 于 2007-6-5 10:42 编辑 ]

风花雪月

程序中没有Sx啊

行者无疆TJ

　由于是离散的,Sx(f)离散成S(k)
IF (F(K+1).GE.8.0D0) THEN
     　　　   IF (F(K+1).LE.1.0D1)  THEN
       　　　 S(K+1)=9.0D-2*U
      　　　  ELSE
     　　   　 S(K+1)=0.0D0
     　　 　  END IF
     　　 　  END IF
      　　  　Y=DSQRT(4.0D0*S(K+1)*DELTAF)
不知是否可以?