哪位有功率谱分析的Fortran程序？发我邮箱840646927@qq.com。 万分感谢。。

Garfield

哪位有功率谱分析的Fortran程序？发我邮箱840646927@qq.com。 万分感谢。。

edwardli

fourier(n,x,a,b,c,s,cta)
　         The discrete Fourier spectrum of one-dimensional series x(n).

求一维序列x(n)的离散Fourier谱分析，s(0:m)离散功率谱，c(0:m)振幅谱，cta(0:m)位相谱，其中m=[n/2.]。

c-----*----------------------------------------------------6---------7--
c     The discrete Fourier spectrum of one-dimensional series x(n).
c     input: n,x(n)
c       n: number of data
c       x(n): raw series
c     output: a(0:m),b(0:m),c(0:m),cta(0:m) and s(0:m)
c       a(0:m): the Fourier coefficient of sine component, m=[n/2.]
c       b(0:m): the Fourier coefficient of cosine component, m=[n/2.]
c       c(0:m): the amplitude spectrum = sqrt(a**2+b**2)/2
c       cta(0:m): the phase angle spectrum =atan(b/a)
c       s(0:m): the discrete power spectrum = (a**2+b**2)/2
c     By Dr. Li Jianping, May 10, 1998.   
      subroutine fourier(n,x,a,b,c,s,cta)
      dimension x(n),a(0:n),b(0:n),c(0:n),cta(0:n),s(0:n)
      pi=3.1415926
      om=2.*pi/n
      m=int(n/2.)
      an=2./n
      do 5 k=0,n
        a(k)=0.
        b(k)=0.
        c(k)=0.
        s(k)=0.
        cta(k)=0.
  5   continue
      call meanvar(n,x,ax,sx,vx)
      a(0)=ax
      b(0)=0.
      do 10 k=1,m
        a(k)=0.
        b(k)=0.
        omg=om*k
        do i=1,n
          i1=i-1
          a(k)=a(k)+x(i)*cos(omg*i1)
          b(k)=b(k)+x(i)*sin(omg*i1)
        enddo
        a(k)=a(k)*an
        b(k)=b(k)*an
  10  continue
      do 20 k=0,m
        sk=a(k)*a(k)+b(k)*b(k)
        c(k)=sqrt(sk)/2.
        s(k)=sk/2.
  20  continue
      cta(0)=0.
      do k=1,m
        cta(k)=atan(b(k)/a(k))
      enddo
      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

edwardli

cspectrum(n,m,x,ol,tl,sl,st95,strw)
　         Subroutine for continuous spectrum analysis of an one-dimensional series x(i) (i=1,...,n).
一维序列x(n)的连续功率谱分析，ol(0:m)频率，tl(0:m)周期，sl(0:m)连续功率谱，st95(0:m)红噪音或白噪音谱的95%置信上限，strw(0:m) 红噪音或白噪音的谱密度，其中m=[n/2.]。


c-----*----------------------------------------------------6---------7--
c     Subroutine for continuous spectrum analysis of an one-dimensional series
c         x(i) (i=1,...,n).
c     Input parameters and arrays: n,m, x(n),
c       n: number of data
c       m: biggest lag time length. Generally, m is between n/3 and n/10.
c       x(n): raw series
c     Output variables:
c        ol(0:m): frequency array
c        tl(0:m): periodic array  
c        sl(0:m): power spectrum. Its sum from 0 to m is 1.
c      st95(0:m): 95% confidence upper limit of red or white noise spectrum
c      strw(0:m): spectrum density of red or white noise
c     By Dr. Li Jianping, May 12,1998.   
      subroutine cspectrum(n,m,x,ol,tl,sl,st95,strw)
      dimension x(n)
      dimension sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
      dimension r(0:m)            ! Work array
      pi=3.1415926
      call autocorrelation(n,m,x,r)
      do 10 l=0,m
        sl(l)=r(0)
        do i=1,m-1
          cc=pi*i/float(m)
          sl(l)=sl(l)+r(i)*(1.+cos(cc))*cos(l*cc)
        enddo
  10  continue
      sl(0)=sl(0)/2.
      sl(m)=sl(m)/2.
      do 20 l=0,m
        sl(l)=sl(l)/m
  20  continue
      do 30 l=1,m
        ol(l)=pi*l/m
        tl(l)=2.*m/l
  30  continue
      ol(0)=0
      tl(0)=100.*tl(1)
      call cnoisetest(n,m,r(1),sl,st95,strw)
      return
      end
c-----*----------------------------------------------------6---------7--
c     This is a subroutine for calculating the lag autocorrelations
c        r(m) of a series x(i).
c     input: n,m,and x(n)
c        n: number of the data x(n)
c        x: the raw series
c        m: maximum lagged length;  m: maximum lag time
c     output: r(0:m)
c        r: lag autocorrelations from 0 to m.
c     By Dr. LI Jianping, April 3, 1998.
      subroutine autocorrelation(n,m,x,r)
      dimension x(n),r(0:m)
      call meanvar(n,x,ax,sx,vx)
      r(0)=1.
      do i=1,n
        x(i)=x(i)-ax
      enddo
      do 10 i=1,m
        sxy=0.
        do j=1,n-i
          sxy=sxy+x(j)*x(i+j)
        enddo
        sxy=sxy/float(n)
        r(i)=sxy/vx
  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     Noise test for continuous spectrum (95% confidence level).
c     strw: spectrum density of red or white noise
c     chi**2 distribution with degrees of freedom v=(2*n-m/2)/m is used.
c     By Dr. Li Jianping, May 12,1998.   
      subroutine cnoisetest(n,m,r1,sl,st95,strw)
      dimension sl(0:m),st95(0:m),strw(0:m)
      dimension chi95(30),chi99(30)
      call chi_table(chi95,chi98,chi99)
      pi=3.1415926
      as=0.
      do 10 l=0,m
        as=as+sl(l)
  10  continue
      as=as/float(m+1)
      r2=r1*r1
      v=(2.*n-m/2.)/float(m)
      if(v.gt.30.)then  !! invaliable from mathmatic table for v>30
        write(*,*)'Unable continue... Please reinput m !'
        stop
      endif
      iv=int(v)
      c95=chi95(iv)+(v-iv)*(chi95(iv+1)-chi95(iv))
      if(r1.gt.0.)then
c       tname='red noise test'
        do l=0,m
          r3=1+r2-2*r1*cos(l*pi/m)
          strw(l)=as*(1.-r2)/r3
          st95(l)=strw(l)*c95/v
        enddo
      else
c       tname='white noise test'
        do l=0,m
          strw(l)=as
          st95(l)=as*c95/v
        enddo
      endif
      return
      end
c-----*----------------------------------------------------6---------7--
c     chi-square distribution
c     Chi-table with right tail probabilities
c       P(X**2>=Xa**2)=a
c       where a (alpha) significance level and (1-a)*100% is confindence level.
c     chi95: chi-square distribution at 95% confidence level.
c     chi98: chi-square distribution at 98% confidence level.
c     chi99: chi-square distribution at 99% confidence level.
c     By Dr. Li Jianping, May 12,1998.   
      subroutine chi_table(chi95,chi98,chi99)
      dimension chi95(30),chi98(30),chi99(30)
      dimension c95(30),c98(30),c99(30) !work array
      data c95/ 3.841, 5.991, 7.815, 9.488,11.070,
     *         12.592,14.067,15.507,16.919,18.307,
     *         19.675,21.026,22.362,23.685,24.996,
     *         26.296,27.587,28.869,30.144,31.410,
     *         32.671,33.924,35.172,36.415,37.652,
     *         38.885,40.113,41.337,42.557,43.773/
      data c98/ 5.412, 7.824, 9.837,11.668,13.388,
     *         15.033,16.622,18.168,19.679,21.161,
     *         22.618,24.054,25.472,26.873,28.259,
     *         29.633,30.995,32.346,33.687,35.020,
     *         36.343,37.659,38.968,40.270,41.566,
     *         42.856,44.140,45.419,46.693,47.662/
      data c99/ 6.635, 9.210,11.345,13.277,15.086,
     *         16.812,18.475,20.090,21.666,23.209,
     *         24.725,26.217,27.688,29.141,30.578,
     *         32.000,33.409,34.805,36.191,37.566,
     *         38.932,40.289,41.638,42.980,44.314,
     *         45.642,46.963,48.278,49.588,50.892/
      do i=1,30
        chi95(i)=c95(i)
        chi98(i)=c98(i)
        chi99(i)=c99(i)
      enddo
      return
      end

wode_q_x

李建平老师的主页上都有，可以下载的。

漏网之鱼

wode_q_x 发表于 2012-4-5 08:57

                               
登录/注册后可看大图

李建平老师的主页上都有，可以下载的。

对头！！李老师的主页超给力！

niuer66

漏网之鱼 发表于 2012-5-5 16:56

                               
登录/注册后可看大图

对头！！李老师的主页超给力！

可以分享下连接么，谢谢了

张_谧

wode_q_x 发表于 2012-4-5 08:57

                               
登录/注册后可看大图

李建平老师的主页上都有，可以下载的。

是啊。李建平老师的主页上都有