星谷胜 三角级数模拟随机序列疑问

行者无疆TJ

最近参考星谷胜的<随机振动分析>一书,3.2.4按三角级数模型的方法(4)中表3-1例3为w在50.24-62.8间 Sx(w)=0.09
         

,而反算出来的谱（图3-15）8-10HZ纵坐标为1.1左右，不知道可不可以理解为：由Sx(f)=Sx(w)*2*pi,再由于表3-1中画的是双边谱，图3-15为单边谱，故再乘2，则可得纵坐标为1.13左右？

然而我自己在反算谱的时候却总是得不到图3-15相近的纵标，我是先通过星谷胜一书3.2.4获取随机序列，再通过一个程序计算谱（计算的是均方根谱），最后得到的谱图纵坐标在3.2左右.

怀疑是随机序列的程序出了错，但花了好长时间也没解决，不知道哪能找现现成的？

行者无疆TJ

没有系统学过随机振动，怀疑是随机序列的程序出了错，但花了几个月时间也没解决，我用的是FORTRAN，我自己所编的获取伪随机序列程序如下，敬请指教！！！！！！！
   也可给我指明另一方法（如MATLAB）不胜感激！！！！！
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

无水1324

不错的东西，没有用过F语言、对随机振动也没有接触过，只能帮你顶一下，让高手帮你解决

行者无疆TJ

有点遗憾，不过还是很感激的！！！！:@)

无水1324

原帖由 行者无疆TJ 于 2007-6-4 13:35 发表
有点遗憾，不过还是很感激的！！！！:@)

是呀
都是学“专业”的所以有些东西肯定不全面

行者无疆TJ

但是我的其实不是太专业化呀,纯粹是随机振动知识.
实质问题就是想把以下有限带宽白噪声的时程序列模拟出来
其中下限频率为8HZ(w为50.24),上限频率是10HZ,Sx(w)=0.09

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