00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 subroutine powerspectrum(t,x,n,iout)
00021 c
00022 c Replace x by its power spectrum, t by frequency.
00023 c Note that it is assumed that t is uniformly spaced.
00024 c
00025 c Since we will usually want to take the logarithm on return,
00026 c SUPPRESS the zero-frequency component.
00027 c
00028 real t(n),x(n)
00029 c
00030 if (n.le.1) return
00031 c
00032 c Discard data to get to a power of 2.
00033 c
00034 i2 = log(float(n))/log(2.)
00035 n = nint(2.**float(i2))
00036 if (iout.ne.0) write(6,*)'Array length truncated to ',n,
00037 & ' for FFT'
00038 c
00039 c Determine the frequency (note storage convention):
00040 c
00041 dti = 1./(t(n) - t(1))
00042 do i=1,n/2
00043 t(i) = i*dti
00044 end do
00045 c
00046 c Take out the mean signal.
00047 c
00048 sum = 0.
00049 do i=1,n
00050 sum = sum + x(i)
00051 end do
00052 c
00053 sum = sum/n
00054 do i = 1,n
00055 x(i) = x(i) - sum
00056 end do
00057 c
00058 c Fourier transform (note packing of result).
00059 c
00060 c Initial data are x(i), i = 1,...,n
00061 c Results are x(1) = F(0), x(2) = F(n/2), x(3) + i x(4) = F(1), etc.
00062 c
00063 call realft(x,n/2,1)
00064 c
00065 c Replace x by the power spectrum.
00066 c
00067 c Zero-frequency component is x(1), note, but suppressed here:
00068 c
00069 c x(1) = x(1)**2
00070 c
00071 x2 = x(2)**2
00072 c
00073 do 50 k=2,n/2
00074 x(k-1) = x(2*k-1)**2 + x(2*k)**2
00075 50 continue
00076 c
00077 x(n/2) = x2
00078 n = n/2
00079 c
00080 99999 end
00081
00082
00083 subroutine realft(data,n,isign)
00084 c
00085 c Routine from Numerical Recipes.
00086 c
00087 real*8 wr,wi,wpr,wpi,wtemp,theta
00088 dimension data(2*n)
00089 c
00090 theta=6.28318530717959d0/2.0d0/dfloat(n)
00091 c1=0.5
00092 if (isign.eq.1) then
00093 c2=-0.5
00094 call four1(data,n,+1)
00095 else
00096 c2=0.5
00097 theta=-theta
00098 endif
00099 wpr=-2.0d0*dsin(0.5d0*theta)**2
00100 wpi=dsin(theta)
00101 wr=1.0d0+wpr
00102 wi=wpi
00103 n2p3=2*n+3
00104 do 95000 i=2,n/2+1
00105 i1=2*i-1
00106 i2=i1+1
00107 i3=n2p3-i2
00108 i4=i3+1
00109 wrs=sngl(wr)
00110 wis=sngl(wi)
00111 h1r=c1*(data(i1)+data(i3))
00112 h1i=c1*(data(i2)-data(i4))
00113 h2r=-c2*(data(i2)+data(i4))
00114 h2i=c2*(data(i1)-data(i3))
00115 data(i1)=h1r+wrs*h2r-wis*h2i
00116 data(i2)=h1i+wrs*h2i+wis*h2r
00117 data(i3)=h1r-wrs*h2r+wis*h2i
00118 data(i4)=-h1i+wrs*h2i+wis*h2r
00119 wtemp=wr
00120 wr=wr*wpr-wi*wpi+wr
00121 wi=wi*wpr+wtemp*wpi+wi
00122 95000 continue
00123 if (isign.eq.1) then
00124 h1r=data(1)
00125 data(1)=h1r+data(2)
00126 data(2)=h1r-data(2)
00127 else
00128 h1r=data(1)
00129 data(1)=c1*(h1r+data(2))
00130 data(2)=c1*(h1r-data(2))
00131 call four1(data,n,-1)
00132 endif
00133 c
00134 end
00135
00136
00137 subroutine four1(data,nn,isign)
00138 c
00139 c Routine from Numerical Recipes.
00140 c
00141 real*8 wr,wi,wpr,wpi,wtemp,theta
00142 dimension data(2*nn)
00143 c
00144 n=2*nn
00145 j=1
00146 do 95000 i=1,n,2
00147 if(j.gt.i)then
00148 tempr=data(j)
00149 tempi=data(j+1)
00150 data(j)=data(i)
00151 data(j+1)=data(i+1)
00152 data(i)=tempr
00153 data(i+1)=tempi
00154 endif
00155 m=n/2
00156 1 if ((m.ge.2).and.(j.gt.m)) then
00157 j=j-m
00158 m=m/2
00159 go to 1
00160 endif
00161 j=j+m
00162 95000 continue
00163 mmax=2
00164 2 if (n.gt.mmax) then
00165 istep=2*mmax
00166 theta=6.28318530717959d0/(isign*mmax)
00167 wpr=-2.d0*dsin(0.5d0*theta)**2
00168 wpi=dsin(theta)
00169 wr=1.d0
00170 wi=0.d0
00171 do 95001 m=1,mmax,2
00172 do 95002 i=m,n,istep
00173 j=i+mmax
00174 tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
00175 tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
00176 data(j)=data(i)-tempr
00177 data(j+1)=data(i+1)-tempi
00178 data(i)=data(i)+tempr
00179 data(i+1)=data(i+1)+tempi
00180 95002 continue
00181 wtemp=wr
00182 wr=wr*wpr-wi*wpi+wr
00183 wi=wi*wpr+wtemp*wpi+wi
00184 95001 continue
00185 mmax=istep
00186 go to 2
00187 endif
00188 c
00189 end
00190
