行者无疆TJ 发表于 2007-6-1 15:52

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

最近参考星谷胜的<随机振动分析>一书,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 发表于 2007-6-3 12:28

星谷胜 三角级数模拟随机序列求助

没有系统学过随机振动，怀疑是随机序列的程序出了错，但花了几个月时间也没解决，我用的是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 50I=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 发表于 2007-6-3 14:32

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

行者无疆TJ 发表于 2007-6-4 13:35

回复 #3 无水1324 的帖子

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

无水1324 发表于 2007-6-4 18:19

本帖最后由 VibInfo 于 2016-5-12 14:50 编辑

原帖由 行者无疆TJ 于 2007-6-4 13:35 发表
有点遗憾，不过还是很感激的！！！！:@)
是呀
都是学“专业”的所以有些东西肯定不全面

行者无疆TJ 发表于 2007-6-5 10:04

回复 #5 无水1324 的帖子

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

[ 本帖最后由 行者无疆TJ 于 2007-6-5 10:20 编辑 ]
页: [1]
查看完整版本: 星谷胜 三角级数模拟随机序列疑问