急！！毕业求助！！关于涨落对耦合细胞中输出信噪比影响

长安云瑾 · 发表于 2014-4-22 11:19:18

我这里需要求200个耦合细胞中每个细胞输出信噪比SNR的平均值，怎么编写程序？同时是否可以在前面构建一个循环求300个，400个、、、1000个的平均值？相关程序如下：

[Fortran] 纯文本查看 复制代码

!! Program to study the signal propagation in Hopf Media ( CSTR This File ) 
Program HBP_ICC_Main
Use MSIMSLMS
Use MSIMSLSS
Use MSFLIB
Use PORTLIB
Use Mymath_Module
Use Graph_Module
Use Model_Module
!! ------------------------------------------------------------------------------- !!
*4
Integer ji_Show( NS ) Integer, Parameter :: NS = 200, NV=3, M_Data = 50000---怎么构建循环，求300个，400个、、、1000个的？不需要多次重复输入。

Real x_Data( NS, M_Data ), x_Data1( M_Data ), x_PSD( 8193 )
Real x_Ary( NS, NV ), x_Ary0( NS, NV ), x_Gauss( NS )
Real xt(NV), dxdt(NV)

Integer ix, ix0, iy(NS), iy0(NS)
Real    xx(NS)
Type (xyCOORD) xyPos
Character a1*3,a2*4 ,a3

Integer, Parameter :: N_C = 10, N_N = 20           !the number of coupling constants and noise
Real  x_Matrix( N_C, N_N, NS, 4 ), xa_tmp( NS )
Real  x_Result(7), xa_R( NS ), xa_Noise( N_N ), xa_CC( N_C )

!! Initialize the Screen for Graphic Drawing ------------------
Call InitGraphWindow( #FF0000 )
Call SetTextWindow(  2, 50, 80, 120 )
jfbk = SetTextColorRGB( 0 )
jfbk = SetColorRGB( #00FFFF )

!Call SetGraphWindow( 10, 10, 650, 510)

! For CSTR, T~2.0; So for DT=1E-3, T~2000*DT ; Set n_Tao=40, we have 50 points
! each period. For NT=2E6, we have 1000 periods. For 40000 data points, the 
! total steps is 40000*40 = 1.6E6. So set n_Try=2E5, NT=2E6. 

! For ICC, T~2.0; for DT=2E-3, n_Tao = 10, we have 100 points per period.
! Set NT=1E6, n_Try = 2E5, n_Data=40000, totally 400 periods ------

DT = 2E-3 ; NT = 6E5 ; n_tao = 10 ; n_try = 1E5 

jf_ShowTime = 0 ; jf_Noise = 1 ; n_Delay = 1;jf_SaveTime = 0           !we can have a result via varing n_Delay
x_CP = .60 ; x_Max = 1.0 ; x_Min = 0.0

jf_BiDirection = 0     ! Birectional Flag        !0(one_way),1(two_way)

n_Show = MIN( NS, 4 ) ; ji_Show(1:n_Show) = (/ 1, 2, 5, 10 /)!, 20, 30 /)

!! ------------------------------------------------------------------------------- !!
n_CC = 1 ; n_Noise = 11
! x_CC = 0.25 ; x_Noise = 0.04 ;        

Do j_C = 1, n_CC     ! For Coupling Constant
Do j_N = 1, n_Noise         ! For Noise Intensity
      !2.4
xa_CC( j_C ) = 3.0 + ( j_C-1 ) *2.0; xa_Noise( j_N ) = -2.0 + (j_N-1) * 0.2  

x_CC = xa_CC( j_C ) ; x_Noise =  10.**( xa_Noise( j_N )        )         ! ~0.004 - 0.63

Call SetGraphWindow( 10, 10, 410, 810 )
j_t1 = Time() ! Start Time For Each Run -----

n_Step = 0 ; x_Gauss( 1:NS ) = 0.0 ;  x_Ary = 0.1 ; n_Data = 0

!!----------------open a file to save the time series--------------------
if( jf_SaveTime == 1  ) then
   write( a1, '(i3.3)' ) NS
   write( a2, '(f4.1)' ) xa_Noise(j_N)
   Write( a3, '(f4.1)' ) x_CC
   Open( 5, file='Data\N_'//a1//'_noise_'//a2//'_CC_'//a3//'.dat', status='unknown' )

endif
!!---------------------------------------------------------
Do it = 1, NT
   if( jf_Noise > 0 .AND. MOD( it-1, n_Delay ) == 0 ) then
       Call RNNOA( jf_Noise, x_Gauss( 1:jf_Noise ) )  
           Do jj = 1, jf_Noise
              if( x_Gauss(jj) < -5.0 ) x_Gauss(jj) = -5.0           ! To Avoid Overflow
                  if( x_Gauss(jj) > 5.0  ) x_Gauss(jj) = 5.0
           Enddo
   endif
   x_Ary0 = x_Ary
   Do js = 1, NS
      ! Using the ICC Model -------------------------------
          xt(1:NV) = x_Ary( js, 1:NV ) 
          xt_cp = x_CP * ( 1.+ x_Noise * x_Gauss(js) )    ! Mulplicative Noise to CP
          Call Derivative( ICC_MODEL, it*DT, xt, dxdt, xt_CP )
      if( js > 1 ) then
              x_Ary(js,1) = x_Ary(js,1) + ( x_Ary0(js-1,1)-x_Ary0(js,1) ) * x_CC * DT   ! Coupling
                  !x_Ary(js,1) = x_Ary(js,1) + x_Ary0(js-1,1 ) * x_CC * DT   ! Coupling
          endif
          if( jf_BiDirection == 1 .AND. js < NS ) then
              x_Ary(js,1) = x_Ary(js,1) + ( x_Ary0(js+1,1)-x_Ary0(js,1) ) * x_CC * DT   ! Coupling
          endif
          !! Now Add the Effect of Noise ----------
          !x_Ary( js,1 ) = x_Ary( js,1 ) + SQRT( 2.*x_Noise*DT ) * x_Gauss(js)
          !x_Ary( js,1 ) = x_Ary( js,1 ) + x_Noise * x_Gauss(js) * DT                          ! Additive Noise 
          x_Ary( js, : ) = x_Ary( js, : ) + dxdt(:) * DT
          Do j_V = 1, NV
             if( x_Ary( js, j_V ) < 0 ) x_Ary( js, j_V ) = 0.0
                 if( x_Ary( js, j_V ) > 1 ) x_Ary( js, j_V ) = 1.0
          Enddo
   Enddo
   if( it <= n_try ) cycle    ! bypass the transient time
   n_Step = n_Step + 1
   if( n_Step == n_tao ) then
           n_data = n_data + 1
           x_Data( 1:NS, n_data ) = x_Ary(1:NS, 1)           
           !! -- Draw the graph on screen ---------------
           if( jf_ShowTime == 1 ) then
                 ix = INT2(REAL(it-n_Try) / (NT-n_Try) * 410.)
                         Do ii = 1, n_Show
                            xx( ii ) = x_Ary( ji_Show(ii), 1 )
                                  iy(ii) = 100*ii - INT2( REAL(xx(ii)-x_Min ) / (x_Max-x_Min) *100. )
                         Enddo                         
                         if( n_Data > 1 ) then
                             Do ii = 1, n_Show
                                    Call Moveto( ix0, iy0(ii), xyPos )
                                    jfbk = LineTo( ix, iy(ii) )
                                 Enddo
                         endif
                         ix0 = ix ; iy0 = iy
                endif
   endif
   if( n_Data >= M_Data ) goto 200
Enddo
200 Continue

! Now calculate the SNR Information and Save the Summary ------------------
Do js = 1, NS

   x_Data1( 1:n_Data ) = x_Data( js, 1:n_Data )-----------主要是这里怎么求200个耦合细胞输出信噪比的平均值！
  ! 20 piece
   Call Get_SNR_Data_2( n_Data, x_Data1, DT*n_Tao, 20, x_Result(1:7) )
   x_Matrix( j_C, j_N, js, 1 ) = x_Result(1)   ! Freq
   x_Matrix( j_C, j_N, js, 2 ) = x_Result(2)   ! Signal
   x_Matrix( j_C, j_N, js, 3 ) = x_Result(3)   ! Noise
   x_Matrix( j_C, j_N, js, 4 ) = x_Result(6)   ! SNR
   xa_R( js ) = x_Result(6)
Enddo