c-----*----------------------------------------------------6---------7--
c     Subroutine for discrete spectrum analysis of an one-dimensional series
c         x(i) (i=1,...,n). 
c     Input parameters and arrays: n,lw,x(n),
c       n: number of data
c       x(n): raw series
c       lw: maximum wave number (=[n/2]).
c     Output variables:
c       ol(lw): frequency array
c       tl(lw): periodic array  
c       sl(lw): power spectrum.
c      st95(lw): 95% confidence upper limit of red or white noise spectrum
c     By Dr. Li Jianping, September 6, 2001.    
      subroutine dspectrum(n,lw,x,tl,sl,st95)
      dimension x(n)
      dimension sl(lw),tl(lw),st95(lw)
      dimension al(lw),bl(lw) !work array
      pi=3.1415926
      c1=2.*pi/float(n)
      do 10 l=1,lw
        al(l)=0.
        bl(l)=0.
        do i=1,n
          al(l)=al(l)+x(i)*cos(c1*l*(i-1))
          bl(l)=bl(l)+x(i)*sin(c1*l*(i-1))
        enddo
        al(l)=2.*al(l)/float(n)
        bl(l)=2.*bl(l)/float(n)
        sl(l)=(al(l)**2+bl(l)**2)/2.
        tl(l)=float(n)/float(l)
 10   continue
      call dnoisetest(n,lw,x,st95)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Significant Test for discrete spectrum (95% confidence level).
c     F distribution test used here, its degrees of freedom numerator (dfn)
c          are 2 and degrees of freedom denominator (dfd) are n-2-1.
c     By Dr. Li Jianping, September 6, 2001.    
      subroutine dnoisetest(n,lw,x,st95)
      dimension x(n),st95(lw)
      dimension fsd95(1000)
      call F_table(fsd95)
      if((n-3).gt.1000)then
        f95=fsd95(1000)
      else
        f95=fsd95(n-3)
      endif
      call meanvar(n,x,ax,sx,vx)
      s=2.*f95*vx/(n-3.+2.*f95)
      do 10 i=1,lw
        st95(i)=s
 10   continue
      return
      end 
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax, standard deviation sx
c       and variance vx of a series x(i) (i=1,...,n).
c     input: n and x(n)
c       n: number of raw series
c       x(n): raw series
c     output: ax, sx and vx
c       ax: the mean value of x(n)
c       sx: the standard deviation of x(n)
c       vx: the variance of x(n)
c     By Dr. LI Jianping, May 6, 1998.
      subroutine meanvar(n,x,ax,sx,vx)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end
c-----*----------------------------------------------------6---------7--
c     The F distribution is a right-skewed distribution used most commonly in 
c         Analysis of Variance. The F distribution is a ration of two Chi-square
c         distributions, and a specific F distribution is denoted by the degrees
c         of freedom for thw numerator Chi-square and the degrees of freedom
c         for the denominator Chi-square.
c         F=S1**2/S2**2,
c         P(F>=Fa)=0.05
c     F(n1,n2): F distribution at 95% confidence level.
c       n1 the numerator degrees of freedom
c       n2 the denominator degrees of freedom
c     Here fsd95(n)=F(2,n) for discrete power spectrum test.
c     By Dr. Li Jianping, September 6, 2001.    
      subroutine F_table(fsd95)
      dimension fsd95(1000)
      dimension fsd(30)
      data fsd/200.00,19.00,9.55,6.94,5.79,5.14,4.74,4.46,4.26,4.10,
     *           3.98, 3.89,3.81,3.74,3.68,3.63,3.59,3.55,3.52,3.49,
     *           3.47, 3.44,3.42,3.40,3.39,3.37,3.35,3.34,3.33,3.32/
      do i=1,30
        fsd95(i)=fsd(i)
      enddo
      do i=31,40
        fsd95(i)=3.32-(3.32-3.23)*(i-30.)/10.
      enddo
      do i=41,50
        fsd95(i)=3.23-(3.23-3.18)*(i-40.)/10.
      enddo
      do i=51,60
        fsd95(i)=3.18-(3.18-3.15)*(i-50.)/10.
      enddo
      do i=61,70
        fsd95(i)=3.15-(3.15-3.13)*(i-60.)/10.
      enddo
      do i=71,80
        fsd95(i)=3.13-(3.13-3.11)*(i-70.)/10.
      enddo
      do i=81,90
        fsd95(i)=3.11-(3.11-3.10)*(i-80.)/10.
      enddo
      do i=91,100
        fsd95(i)=3.10-(3.10-3.09)*(i-90.)/10.
      enddo
      do i=101,500
        fsd95(i)=3.09-(3.09-3.01)*(i-100.)/400.
      enddo
      do i=501,1000
        fsd95(i)=3.01-(3.01-3.00)*(i-500.)/500.
      enddo
      return
      end