program mtmsvdrecon c c Copyright (C) 1995 c c Michael Mann c c perform spatial & temporal mode reconstructions at a particular c frequency for a peak detected by MTM-SVD method c c v1.0 c c uses stripped version of original multitaper spectral code of c Jeffrey Park (c) c c mmax - maximum number of spatial dofs desired c (ie, max # of time series) c ntapmax - maximum number of spectral dofs desired c nscanmax - maximum length of time series c nitermax - maximum number of independent windows c calculated in evolutive analysis c parameter (mmax=25, $ nscanmax=1200,ntapmax=3,nitermax=400) character*15 filename1,filename2,filename3 character *19 filename4 character*1 cfnum,cf1,cf2,cf3,cf4 integer igrid(mmax) real aaa(nscanmax,mmax),anew(nscanmax,mmax) real alow(nscanmax,mmax) real yin(nscanmax),yout(nscanmax) real rcon(mmax,nscanmax),recon_project(mmax,nscanmax) real rconbest(mmax,nscanmax) real recon_proj1(mmax,nscanmax), $ recon_proj2(mmax,nscanmax),recon_proj3(mmax,nscanmax) real rcon1(mmax,nscanmax), $ rcon2(mmax,nscanmax),rcon3(mmax,nscanmax) integer ncounts(nscanmax) real cmag(mmax),angle(mmax),sum1(mmax), $ sum2(mmax) real *8 ar,ai,S(ntapmax),Save(ntapmax), $ partvar(ntapmax) real *8 envmax(ntapmax),env1(2,8194),env2(2,8194), $ env3(2,8194) real envmag1,envmag2,envmag3 real*8 ell,env0,cs,sn,snprm real dum,sd(mmax),weight(mmax),datave(mmax) real snorm(mmax) c integer ifreq(nscanmax) real afreq(nscanmax),tasold(8192,ntapmax) real lambda1,lambda2,lambda3, $ lambmin1,lambmin2,lambmin3,sum,amp real *8 aveamp(mmax,mmax), avephase(mmax,mmax) complex *16 AA(mmax,ntapmax),U(mmax,mmax), $ V(ntapmax,ntapmax) complex *16 Uave(mmax,mmax),Vave(ntapmax,ntapmax) c c character*4 chead(158) complex zed c common/npiprol/anpi common/nnaamm/iwflag common/work/b1(8194) common/taperz/ta(4100,16),tai(4100,16),dcf(4100,16),amu(4100,4) common/staple/tad(8200) common/stap2/tas(8192,16) common/stap4/cs(8200),sn(8200) common/envel1/ell(20) common/data/a(50000) c common/pltopt/npad,nst,mpts,nmin,nmax,t(3),fmin,fmax,nf, x ddt,nwin,inc,tt(3) c equivalence (zed,reim),(iy,dum),(iah,ahead),(iah,chead), x (tad,env0) c c dimension dum(35),ahead(158),env0(2,1000) iwflag=0 zero = 0.0 one = 1.0 two = 2.0 half = 0.5 pi=3.14159265358979 c c c dt = time interval spacing in units of years c nscan = length of time series in units of "dt" c iabv = number of time series in analysis (ie, "spatial" array) c c c c dt = time interval spacing in units of years c nscan = length of time series in units of "dt" c iabv = number of time series in analysis (ie, "spatial" array) c c presently set for 25 records of 1 month resolution and 100 year c length (** make sure array dimensions are large enough: c "ifmax" must be greater than nscan **) c so nscan = 100*12 = 1200, iabv = 25 c c dt = 1/12 so that time units/frequencies will be in terms c of annual time units c c **** THESE SETTINGS ARE DATASET DEPENDENT!!! **** c dt = 0.0833333 nscan = 1200 iabv = 25 nscan = 1200 c c define initial time as "year 0" c nbegin = 0 c c c initialize uniform weights on dataseries -- the user c may decide to set e.g., areal (ie, cos(latitude) ) weights c on the constituent dataseries, etc. c do i=1,mmax weight(i)=1.0 end do c c c set frequency range for DFT, frequency spacing, zero c padding, etc. c ** NOTE: maximum padded DFT length is npad=8192 ** c npad must be a power of 2 greater than or equal to c nscan. Larger values provide a more interpolated c frequency axis c npad = 4096 c c f1 = minimum frequency for analysis = 0 c f2 = maximum frequency (in cyc/year -- this can be varied up to c the nyquist sampling frequency = 0.5/dt) c c ( presently set for the range 0 to 0.5 cycle year, or, periods greater c than tau=2 years. This is a smaller fraction of the full Nyquist c interval f=0 to f=6.0 cycle/year for monthly sampling) c f1 = 0.0 f2 = 0.5 c c ddf=1./(npad*dt) eps=.1*ddf nf1=(f1+eps)/ddf nf2=(f2+eps)/ddf nf=nf2-nf1+1 nf = nf-1 nfmax = nf c t(1)=f1 t(2)=ddf t(3)=ddf fmin=nf1*ddf fmax=nf2*ddf nmin=nf1*2+1 nmax=nf2*2+2 c c set time interval default parameters c n1 = 1 n2 = nscan neval = n2-n1+1 c c read in data. c c * THIS SECTION WILL HAVE TO MODIFIED DEPENDING ON THE DATA c ANALYZED * c open (unit=11,file='synth1.dat', $ status='old') c open (unit=12,file='synth2.dat', $ status='old') c open (unit=13,file='synth3.dat', $ status='old') c open (unit=14,file='synth4.dat', $ status='old') c open (unit=15,file='synth5.dat', $ status='old') c do j=1,5 k1 = (j-1)*5+1 k2 = k1+4 do i=1,nscan read (10+j,777) adum,(aaa(i,k),k=k1,k2) end do close (unit=10+j) end do 777 format (6f9.3) c c write (6,*) 'read in: ' write (6,*) iabv, ' time series' write (6,*) nscan*dt, ' years of data...' c c c for each time series, subtract mean, normalize by standard c deviation for each month, c also set weights for the different dataseries. c ( for the the present case, all dataseries are equally weighted for c analysis) c c open (unit=15,file='recon-stats.out', $ status='unknown') do j=1,iabv sum=0.0 do i=1,nscan sum=sum+aaa(i,j) end do datave(j)=sum/float(nscan) var = 0.0 do i=1,nscan anew(i,j)=aaa(i,j)-datave(j) alow(i,j)=anew(i,j) var=var + anew(i,j)**2 end do sd(j)=sqrt(var)/float(nscan) snorm(j)=sd(j) do i=1,nscan anew(i,j)=weight(j)*anew(i,j)/sd(j) end do write (15,*) j,snorm(j),datave(j) end do close (unit=15) c c default settings c ievolute=0 iif = 1 ff0 = 0.0 freq=0.0 do i=1,nscan afreq(i) = 0.0 ifreq(i) = 1 end do npi = 2 anpi=float(npi) nwin=3 imode = 0 nmode = 1 iseg = 0 iformat = 0 iref = 0 jsmoo = 0 iaverage = 0 ilambda = 0 lambda1 = 1.0 lambda2 = 0.0 lambda3 = 0.0 nloop1 = 0 nloop2 = 0 irecon = 0 imoderecon = 0 ioutput = 0 ilow = 0 nmove = n2-n1+1 neval = n2-n1+1 npts = neval ninterv = nmove nsites = 0 c 444 write (6,*) 'present settings' write (6,*) '1) use',nwin, ' x ',npi,' pi tapers' if (imode.eq.0) then write (6,*) '2) determine primary mode only' else if (imode.eq.1) then write (6,*) '2) determine primary and secondary mode' else if (imode.eq.2) then write (6,*) '2) determine modes ',1,' through ',nwin-1 endif endif endif if (iseg.le.1) then write (6,*) '3) fixed freq recon: segment ',n1,' to ',n2 else write (6,*) '3) evolutive reconstruction from ', $ n1,' to ',n2 write (6,*) ' window = ',nmove write (6,*) ' interval = ',ninterv endif if (iformat.eq.0) then write (6,*) '4) complex spatial pattern format' else write (6,*) '4) successive real-valued spatial snapshots' if (icycle.eq.0) then write (6,*) ' -- reconstruct 1/2 cycle' else write (6,*) ' -- reconstruct full cycle' endif write (6,*) ' -- ',nsnap,' snapshots' endif write (6,*) '5) reference site for zero/initial phase: #',iref if (jsmoo.eq.0) then write (6,*) '6) minimum norm time-domain inversion' else if (jsmoo.eq.1) then write (6,*) '6) minimum slope time-domain inversion' else if (jsmoo.eq.2) then write (6,*) '6) minimum roughness time-domain inversion' else write (6,*) '6) minimum error time-domain inversion' if (ilambda.ne.1) then write (6,*) ' min norm weight:',lambda1 write (6,*) ' min rough weight:',lambda2 write (6,*) ' min slope wieght:',lambda3 endif endif endif endif c if (iaverage.eq.0) then c write (6,*) '7) do not calculate spatial pattern averages' c else c write (6,*) '7) calculate spatial pattern averages' c endif if (irecon.eq.0) then write (6,*) '7) no specific time-reconstructions' else write (6,*) '7) perform time reconstructions' write (6,*) ' mode #',imoderecon write (6,*) ' sites #:' if (nsites.eq.iabv) then write (6,*) ' all ',iabv,' sites' else write (6,1010) (igrid(i),i=1,nsites) endif endif if (ievolute.eq.0) then write (6,*) '8) fixed-frequency reconstruction' write (6,*) ' frequency = ',freq else write (6,*) '8) variable frequency reconstruction' write (6,*) " read data from file 'frequency.dat'" endif 1010 format (3x,25i3) c write (6,*) '-----------------' write (6,*) 'item to change (0=continue)' read (5,*) icont if (icont.eq.1) then goto 1111 else if (icont.eq.2) then goto 2222 else if (icont.eq.3) then goto 3333 else if (icont.eq.4) then goto 4444 else if (icont.eq.5) then goto 5555 else if (icont.eq.6) then goto 6666 else c if (icont.eq.7) then c goto 7777 c else if (icont.eq.7) then goto 8888 else if (icont.eq.8) then goto 9999 else goto 9876 endif endif c endif endif endif endif endif endif endif c c anpi = desired time/frequency bandwith "p" (ie, "p pi tapers") c nwin = number of tapers used (maximum is 2p-1) c 1111 write (6,*) 'time-frequency bandwidth product "p"' read (5,*) npinew if (ntapmax.lt.2*npinew-1) then write (6,*) 'sorry--exceeds maximum bandwidth' else write (6,*) 'number of tapers' read (5,*) nwinnew if (nwinnew.gt.2*npinew-1) then write (6,*) 'sorry--cannot exceed 2*p-1' else npi = npinew anpi = float(npi) nwin=nwinnew endif endif goto 444 c c 2222 write (6,*) 'reconstruct' write (6,*) '(0) primary mode' write (6,*) '(1) primary and secondary mode' write (6,*) '(2) modes ',1,' through ',nwin-1 read (5,*) imode nmode=imode+1 if (nmode.gt.2) nmode=nwin-1 goto 444 c c signal can be reconstructed using the c SVD of the entire series, or a particular segment c of the time series (e.g., a statistically signficant c time segment indicated by the evolutive SVD fractional c variance spectrum calculated earlier). Time reconstructions c can be performed EVOLUTIVELY as well...The spatial pattern c returned will be the average over all windows... c 3333 write (6,*) 'reconstruct signal based on' write (6,*) '(0) entire series' write (6,*) '(1) particular time segment' write (6,*) '(2) evolutive reconstruction' read (5,*) iseg if (iseg.eq.1) then write (6,*) 'subsegment between 1 and ',nscan,' :' write (6,*) 'initial time index (eg month or year #)' read (5,*) n1 write (6,*) 'final time index (eg month or year #)' read (5,*) n2 npts = n2-n1+1 nmove = npts ninterv = nmove else if (iseg.eq.0) then n1 = 1 n2 = nscan npts = n2-n1+1 else write (6,*) 'subsegment between 1 and ',nscan,' :' write (6,*) 'initial time index (eg month or year #)' read (5,*) n1 write (6,*) 'final time index (eg month or year #)' read (5,*) n2 write (6,*) 'window width (in # of time samples)' read (5,*) nmove write (6,*) 'evaluation interval (in # of time samples)' read (5,*) ninterv npts = nmove endif endif goto 444 c 4444 write (6,*) 'format for spatial pattern information' write (6,*) '(0) complex pattern (amplitude and phase)' write (6,*) '(1) successive real-valued snapshots' read (5,*) iformat if (iformat.eq.1) then write (6,*) '(0) reconstruct half cycle' write (6,*) '(1) reconstruct full cycle' read (5,*) icycle write (6,*) 'number of snapshots' read (5,*) nsnap endif goto 444 c c for definiteness in pattern reconstruction, zero phase of cycle c needs to be defined by maximum value at particular reference site c 5555 if (iformat.eq.0) then write (6,*) 'need reference site for defining zero phase...' else write (6,*) 'need reference site to define initial snapshot...' endif write (6,*) 'reference point (',1,' to ',iabv,' ) (0=absolute)' read (5,*) iref goto 444 c c "minimum norm" : minimize the overall size of the envelope of the c oscillation. Generally the most faithful constraint, but c may artificially minimize amplitude ends of time interval. c "minimum slope" :minimize the average slope of the envelope of c the oscillation. Particularly useful for reconstructing c secular trends in the series which do not have the periodic c property of non-secular oscillations. c c "minimum roughness" : minimize the average 2nd derivative of the c envelope of the oscillation. Useful if the amplitude of the c envelope is believed to be changing abruptly near the ends of c the time interval c c "minimum error" : minimize the MSE between raw data and reconstruction c over all possible linear combinations of the alternative c a priori constraints described above c 6666 write (6,*) 'inversion for time domain' write (6,*) '(0) minimum norm' write (6,*) '(1) minimum slope' write (6,*) '(2) minimum roughness' write (6,*) '(3) minimum MSE combination' read (5,*) jsmoo if (jsmoo.eq.3) then lambda = 0.0 write (6,*) '(0) specify lambda' write (6,*) '(1) determine lambda by minimization' read (5,*) ilambda 6667 if (ilambda.eq.0) then write (6,*) 'weight on min norm solution (0-1)' read (5,*) lambda1 write (6,*) 'weight on min rough solution (0-1)' read (5,*) lambda2 write (6,*) 'weight on min slope solution (0-1)' read (5,*) lambda3 if ((lambda1+lambda2+lambda3.gt.1.05).or. $ (lambda1+lambda2+lambda3.lt.0.95)) then write (6,*) 'weights must sum to one!!' goto 6667 endif lambda3 = one-lambda1-lambda2 else write (6,*) 'minimize all combinations of' write (6,*) '(0) min norm/min roughness' write (6,*) '(1) min norm/min roughness/min slope' read (5,*) lmin nloop1 = 50 nloop2 = 50*lmin endif endif goto 444 c c if desired, averages of area-weighted c spatial-patterns can be calculated using coordinates c of locations of time series c c7777 write (6,*) 'calculate area-weighted spatial averages' c write (6,*) '(0) No' c write (6,*) '(1) Yes' c read (5,*) iaverage c goto 444 c 8888 write (6,*) 'Do specific site time-reconstructions?' write (6,*) '(0) No' write (6,*) '(1) Yes' read (5,*) irecon if (irecon.eq.1) then write (6,*) 'mode to reconstruct (',1,' to ',nwin-1,')' read (5,*) imoderecon if (imoderecon.gt.nwin-1) imoderecon=nwin-1 if (imoderecon.gt.1) then imode=2 nmode=nwin-1 endif write (6,*) 'reconstruct' write (6,*) '(0) all ',iabv,' sites' write (6,*) '(1) particular site' write (6,*) '(2) particular subset of sites' read (5,*) ireconsite nsites=0 if (ireconsite.eq.1) then nsites=1 write (6,*) 'site # to reconstruct (',1,' to ',iabv,')' read (5,*) igrid(nsites) else if (ireconsite.eq.2) then write (6,*) 'site # to reconstruct (',1, $ ' to ',iabv,')' write (6,*) 'enter "0" to terminate list...' 111 nsites=nsites+1 read (5,*) igrid(nsites) if (igrid(nsites).ne.0) goto 111 nsites = nsites-1 else nsites=iabv do i=1,nsites igrid(i)=i end do endif endif write (6,*) 'output time reconstructions?' write (6,*) '(0) No' write (6,*) '(1) Yes' read (5,*) ioutput write (6,*) 'reference time series' write (6,*) '(0) no lowpass' write (6,*) '(1) interseasonal (6 month) lowpass' write (6,*) '(2) interannual (2 yr) lowpass' read (5,*) ilow fcut = 0.0 if (ilow.eq.1) then fcut=2.0 else if (ilow.eq.2) then fcut = 0.5 endif endif endif goto 444 c c 9999 write (6,*) 'Perform time-reconstruction based on' write (6,*) '(0) fixed frequency' write (6,*) '(1) evolving frequency pattern' read (5,*) ievolute if (ievolute.eq.0) then c c get frequency for SVD, determine associated DFT index c write (6,*) 'frequency (cycles/year)' read (5,*) freq iif = 1+(freq+eps)/ddf freq = ddf*(iif-1) do i=1,nscan ifreq(i)=iif afreq(i)=freq end do else open (unit=4,file='frequency.dat',status='old') read (4,*,end=8765) itold,afreqold ifreq(1)=1+(afreqold+eps)/ddf afreq(1)=ddf*(ifreq(1)-1) do i=1,nscan read (4,*,end=8765) it,afreqnew slope = (afreqnew-afreqold)/float(it-itold) do j=itold+1,it afreq(j)=afreqold+float(j-itold)*slope ifreq(j)=1+(afreq(j)+eps)/ddf afreq(j)=ddf*(ifreq(j)-1) end do afreqold = afreqnew itold = it end do 8765 close (unit=4) open (unit=14,file='frequency.out',status='unknown') freqsum = 0.0 freqmin = 999999.0 freqmax=-freqmin do i=1,nscan write (14,*) i,ifreq(i),afreq(i) if (afreq(i).gt.freqmax) freqmax=afreq(i) if (afreq(i).lt.freqmin) freqmin=afreq(i) freqsum = freqsum + afreq(i) end do freq = freqsum/float(nscan) close (unit=14) write (6,*) 'min freq: ',freqmin write (6,*) 'max freq: ',freqmax write (6,*) 'ave freq: ',freq endif goto 444 c c c 9876 niter = 1 if (iseg.eq.2) then niter = (neval-nmove)/ninterv + 1 if (niter.gt.nitermax) then write (6,*) 'too many iterations!' write (6,*) ' -- decrease step/increase window width' goto 444 endif endif c c c OK...we're happy with the settings. Begin calculations... c c First, output some appropriate information c write (6,*) 'nscan ',nscan write (6,*) 'neval ',neval write (6,*) 'dt ',dt write (6,*) 'nmove ',nmove write (6,*) 'ninterv ',ninterv write (6,*) 'nmode ',nmode c open (unit=9,file='recon-data.out',status='unknown') c c c calculate reference lowpass of each series if requested c c if (ilow.gt.0) then write (6,*) 'determining lowpass filtered series...' do jjj=1,nsites iii=igrid(jjj) do i=1,nscan yin(i)=alow(i,iii) end do call lowpass(yin,yout,nscan,fcut,dt) do i=1,nscan alow(i,iii)=yout(i) end do end do write (6,*) 'done calculating lowpasses...' endif c c c construct tapers appropriate for time interval/window width c call taper(npts,nwin,ell) c c c decimate the tapers c do i=1,npts do k=1,nwin tasold(i,k)=tas(i,k) end do end do c inc=(npts-1)/500+1 finc=float(inc) j=0 do i=1,npts,inc j=j+1 do k=1,nwin tas(j,k)=tas(i,k) end do end do mpts=j ddt=dt*inc c c initialize running averages c do i=1,iabv do j=1,nmode Uave(i,j)=(0.d0,0.d0) aveamp(i,j) = 0.d0 avephase(i,j) = 0.d0 end do do j=1,nscan recon_project(i,j)=zero end do end do do i=1,nwin do j=1,nmode Vave(i,j)=(0.d0,0.d0) end do Save(i)=0.d0 end do do i=1,nscan ncounts(i)=0 end do do i=1,nmode envmax(i) = -999999.0 end do c c set parameters for SVD c M = iabv N = nwin NV = N NU = N IP = 0 c c initialize variance sum for duration of series c sumtot1 = 0.0 sumtot2 = 0.0 summintot1 = 0.0 summintot2 = 0.0 c c c EVOLUTIVE RECONSTRUCTION CASE (reduces to single loop c for non-evolutive case) c if (iseg.eq.2) $ write (6,*) 'beginning 1st of ',niter, ' iteration/s...' c c 10 do it = 1,niter c c nn1 = n1+(it-1)*ninterv nn2 = nn1+npts-1 if ((it.eq.niter).and.(nn2.ne.n2)) then nn2 = n2 nn1 = n2-npts+1 endif ncent = (nn1+nn2-1)/2 if (iseg.eq.2) then open (unit=28,file='recon-status.out',status='unknown') write (28,*) 'it n1 n2 ncent ',it,nn1,nn2,ncent close (unit=28) write (6,*) 'it n1 n2 ncent ',it,nn1,nn2,ncent endif ff0 = afreq(ncent) iif = ifreq(ncent) iwhere = 2*(nf1+iif)-1 c c calculate DFTs for each time series c do ii=1,iabv c demean gridpoint series locally sum = 0. do i=nn1,nn2 sum = sum+anew(i,ii) end do sum = sum/float(nn2-nn1+1) do iwin=1,nwin do i=1,npad b1(i)=0. end do do i=nn1,nn2 b1(i-nn1+1)=(anew(i,ii)-sum)*tasold(i+1-nn1,iwin) c if ((ii.eq.1).and.(iwin.eq.1)) then c write (2,*) i,anew(i,ii)-sum,tasold(i+1-nn1,iwin) c endif end do call realft(b1,npad/2,1) b1(npad-1)=b1(2) b1(npad)=0. b1(2)=0. ar=b1(iwhere) ai=b1(iwhere+1) AA(ii,iwin)=dcmplx(ar,ai) end do end do c c perform SVD for given segment c call CSVD(AA,mmax,nwin,M,N,IP,NU,NV,S,U,V) c c keep running average of spatial, spectral eofs c c if not error minimization option c first, calculate reference phase for given time window c do i=1,iabv do j=1,nmode angle_ref = 180.0*atan2(dimag(U(iref,j)), $ dreal(U(iref,j)))/pi Uave(i,j)=Uave(i,j)+U(i,j) amp = dreal(U(i,j))**2+dimag(U(i,j))**2 aveamp(i,j)=aveamp(i,j)+amp if (amp.gt.zero) avephase(i,j)=avephase(i,j)+ $ (180.0*atan2(dimag(U(i,j)),dreal(U(i,j)))/pi-angle_ref) end do end do do i=1,nwin do j=1,nmode Vave(i,j)=Vave(i,j)+V(i,j) end do Save(i)=Save(i)+S(i) end do c c do i=1,2*nscanmax a(i)=0.0 end do c c MODE LOOP c 123 do imode=1,nmode c c load in spectral envelopes for each mode c do iwin=1,nwin do j=1,nscanmax ta(j,iwin)=0.0 tai(j,iwin)=0.0 end do ta(iif,iwin)=dreal(V(iwin,imode)) tai(iif,iwin)=-dimag(V(iwin,imode)) if (iif.eq.1) ta(iif,iwin)=ta(iif,iwin)/2.0 end do c c invert spectral envelope to get time envelope c call envel(ff0,0) do i=1,mpts env1(1,i)=env0(1,i) env1(2,i)=env0(2,i) end do call envel(ff0,2) do i=1,mpts env2(1,i)=env0(1,i) env2(2,i)=env0(2,i) end do call envel(ff0,1) do i=1,mpts env3(1,i)=env0(1,i) env3(2,i)=env0(2,i) end do call cossin(npts,cs,sn,1.d0,0.d0,dble(ff0),dble(dt)) c c c do time reconstructions for indicated mode for c all sites indicated c c do i=1,mpts-1 ii=(i-1)*inc do j=1,inc c c determine time envelope A_k(t) of mode, and c find maximum amplitude attained over record c im = ii+j a11=(env1(1,i)*(inc+1-j)+env1(1,i+1)*(j-1))/finc a21=(env1(2,i)*(inc+1-j)+env1(2,i+1)*(j-1))/finc a12=(env2(1,i)*(inc+1-j)+env2(1,i+1)*(j-1))/finc a22=(env2(2,i)*(inc+1-j)+env2(2,i+1)*(j-1))/finc a13=(env3(1,i)*(inc+1-j)+env3(1,i+1)*(j-1))/finc a23=(env3(2,i)*(inc+1-j)+env3(2,i+1)*(j-1))/finc envmag1 = sqrt(a11**2+a21**2) envmag2 = sqrt(a12**2+a22**2) envmag3 = sqrt(a13**2+a23**2) if (envmag1.gt.envmax(imode)) envmax(imode)=envmag1 if (envmag2.gt.envmax(imode)) envmax(imode)=envmag2 if (envmag3.gt.envmax(imode)) envmax(imode)=envmag3 c c do specific requested time reconstructions c if (imode.eq.imoderecon) then c c take complex conjugate of carrier frequency sinusoid c snprm = -sn(im) c c take Re{A_k(t) U_l(x) exp(iw_p t)} to reconstruct c oscillation at grid point k c do jjj=1,nsites iii=igrid(jjj) a(npts+im) = $ dreal(U(iii,imode))*(a11*cs(im)-a21*snprm) $ -dimag(U(iii,imode))*(a11*snprm+a21*cs(im)) rcon1(iii,im)=S(imode)*snorm(iii) $ *a(npts+im)/weight(iii) c a(npts+im) = $ dreal(U(iii,imode))*(a12*cs(im)-a22*snprm) $ -dimag(U(iii,imode))*(a12*snprm+a22*cs(im)) rcon2(iii,im)=S(imode)*snorm(iii) $ *a(npts+im)/weight(iii) c a(npts+im) = $ dreal(U(iii,imode))*(a13*cs(im)-a23*snprm) $ -dimag(U(iii,imode))*(a13*snprm+a23*cs(im)) rcon3(iii,im)=S(imode)*snorm(iii) $ *a(npts+im)/weight(iii) end do endif c c done w/ specific time reconstructions c end do end do c c linear extrapolate reconstructions one point c ii=(mpts-1)*inc e11=2.*env1(1,mpts)-env1(1,mpts-1) e21=2.*env1(2,mpts)-env1(2,mpts-1) e12=2.*env2(1,mpts)-env2(1,mpts-1) e22=2.*env2(2,mpts)-env2(2,mpts-1) e13=2.*env3(1,mpts)-env3(1,mpts-1) e23=2.*env3(2,mpts)-env3(2,mpts-1) do j=1,inc im = ii+j a11=(env1(1,mpts)*(inc+1-j)+e11*(j-1))/finc a21=(env1(2,mpts)*(inc+1-j)+e21*(j-1))/finc a12=(env2(1,mpts)*(inc+1-j)+e12*(j-1))/finc a22=(env2(2,mpts)*(inc+1-j)+e22*(j-1))/finc a13=(env3(1,mpts)*(inc+1-j)+e13*(j-1))/finc a23=(env3(2,mpts)*(inc+1-j)+e23*(j-1))/finc c c do specific requested time reconstructions c if (imode.eq.imoderecon) then c c take complex conjugate of carrier frequency sinusoid c snprm = -sn(im) do jjj=1,nsites iii=igrid(jjj) a(npts+im) = $ dreal(U(iii,imode))*(a11*cs(im)-a21*snprm) $ - dimag(U(iii,imode))*(a11*snprm+a21*cs(im)) rcon1(iii,im)=S(imode)*snorm(iii)*a(npts+im) $ /weight(iii) c a(npts+im) = $ dreal(U(iii,imode))*(a12*cs(im)-a22*snprm) $ - dimag(U(iii,imode))*(a12*snprm+a22*cs(im)) rcon2(iii,im)=S(imode)*snorm(iii)*a(npts+im) $ /weight(iii) c a(npts+im) = $ dreal(U(iii,imode))*(a13*cs(im)-a23*snprm) $ - dimag(U(iii,imode))*(a13*snprm+a23*cs(im)) rcon3(iii,im)=S(imode)*snorm(iii)*a(npts+im) $ /weight(iii) end do endif c c done w/ specific time reconstructions c end do c c c END MODE LOOP c end do c c c sum contributions to overall time projection form c reconstructions in overlapping windows c do imloc=1,npts imglobe = imloc+nn1-1 ncounts(imglobe)=ncounts(imglobe)+1 do iii=1,iabv recon_proj1(iii,imglobe) = recon_proj1(iii,imglobe) $ + rcon1(iii,imloc) recon_proj2(iii,imglobe) = recon_proj2(iii,imglobe) $ + rcon2(iii,imloc) recon_proj3(iii,imglobe) = recon_proj3(iii,imglobe) $ + rcon3(iii,imloc) end do end do c c END EVOLUTIVE LOOP c 100 end do if (iseg.eq.2) write (6,*) 'finished time loop...' c c calculate projection filtered reconstructions c c ERROR MINIMIZATION LOOP c summin1 = 1.0e+36 lambmin1 = lambda1 lambmin2 = lambda2 lambmin3 = lambda3 do irep1=1,nloop1+1 do irep2=1,nloop2+1 c if (jsmoo.eq.3) then if (ilambda.eq.1) then lambda1 = float(irep1-1)/float(nloop1) lambda3 = 0.0 if (lmin.gt.0) lambda3 = float(irep2-1)/float(nloop2) lambda2=one-lambda1-lambda3 if (lambda2.lt.0.0) goto 3456 endif else lambda1 = 0.0 lambda2 = 0.0 lambda3 = 0.0 if (jsmoo.eq.0) then lambda1=one else if (jsmoo.eq.1) then lambda3=one else lambda2=one endif endif endif c c determine fractional explained variance as function c of lambda over indicated sites c if (irecon.gt.0) then sumtot1 = 0.0 sumtot2 = 0.0 do jjj=1,nsites iii=igrid(jjj) do i=n1,n2 rcon(iii,i)=lambda1*recon_proj1(iii,i)/float(ncounts(i)) $ + lambda2*recon_proj2(iii,i)/float(ncounts(i)) $ + lambda3*recon_proj3(iii,i)/float(ncounts(i)) term1 = ((rcon(iii,i)-alow(i,iii)) $ *weight(iii)/snorm(iii))**2 term2 = (alow(i,iii)*weight(iii) $ /snorm(iii))**2 sumtot1 = sumtot1+term1 sumtot2 = sumtot2+term2 end do end do c write (6,*) 'lambda1 lambda2 lambda3 ', $ lambda1,lambda2,lambda3 write (9,*) 'lambda1 lambda2 lambda3 ', $ lambda1,lambda2,lambda3 if (ilow.gt.0) then write (6,*) '% of lowpassed variance over all ', $ nsites,' sites: ',100.0*(one-sumtot1/sumtot2) write (9,*) '% of lowpassed variance over all ', $ nsites,' sites: ',100.0*(one-sumtot1/sumtot2) else write (6,*) '% of raw variance over all ', $ nsites,' sites: ',100.0*(one-sumtot1/sumtot2) write (9,*) '% of raw variance over all ', $ nsites,' sites: ',100.0*(one-sumtot1/sumtot2) endif if (sumtot1.lt.summin1) then summin1 = sumtot1 summin2 = sumtot2 summintot1 = sumtot1 summintot2 = sumtot2 lambmin1 = lambda1 lambmin2 = lambda2 lambmin3 = lambda3 do jjj=1,nsites iii=igrid(jjj) do i=n1,n2 rconbest(iii,i)=rcon(iii,i) end do end do endif endif c c c END ERROR MINIMIZATION loop c c 3456 end do end do c c c c open output files c c recon-spat.eof[i] (i=1,...,nmode) : spatial eof_i c recon-spec.eof[i] " : spectral eof_i c recon-spat.pat[i] " : spatial pattern_i c do i=1,nmode cfnum = char(48+i) filename1 = 'recon-spat.eof'//cfnum filename2 = 'recon-spec.eof'//cfnum filename3 = 'recon-spat.pat'//cfnum open (unit=13+i-1,file=filename1, $ status='unknown') open (unit=23+i-1,file=filename2, $ status='unknown') open (unit=33+i-1,file=filename3, $ status='unknown') end do c c determine average spatial/spectral EOFs c do i=1,iabv do j=1,nmode U(i,j)=Uave(i,j)/float(niter) aveamp(i,j)=dsqrt(aveamp(i,j)/float(niter)) avephase(i,j)=avephase(i,j)/float(niter) end do end do do i=1,nwin do j=1,nmode V(i,j)=Vave(i,j)/float(niter) end do S(i)=Save(i)/float(niter) end do c c write out spatial eofs, spectral eofs c do i=1,iabv do j=1,nmode write (13+j-1,*) dreal(U(i,j)),dimag(U(i,j)) end do end do do i=1,nwin do j=1,nmode write (23+j-1,*) dreal(V(i,j)),dimag(V(i,j)) end do end do c c print out any time reconstructions if indicated c c format is: column 1 column2 column3 c time reconstructed mode raw time series c if (ioutput.eq.1) then do jjj=1,nsites iii=igrid(jjj) idig1 = iii/1000 idig2 = (iii-1000*idig1)/100 idig3 = (iii-1000*idig1-100*idig2)/10 idig4 = iii-1000*idig1-100*idig2-10*idig3 cf1 = char(48+idig1) cf2 = char(48+idig2) cf3 = char(48+idig3) cf4 = char(48+idig4) filename4 = 'recon-time.site'//cf1//cf2//cf3//cf4 open (unit=19,file=filename4,status='unknown') do i=n1,n2 write (19,222) float(nbegin)+i*dt,rcon(iii,i), $ alow(i,iii),anew(i,iii)*sd(iii)/weight(iii) end do close (unit=19) end do endif 222 format (f10.2,3f14.4) c c determine spatial patterns (corresponding to peak amplitude c oscillation over duration of record--note this is 1/2 peak-to-peak) c c amplitude_i = Singular value(j) * max|time envelope| c * value of spatial eof(j)_i * standard deviation_i c / weight_i c c |amplitude(i)| = peak amplitude variation at site i c phase(amplitude(i)) = temporal phase (relative to maximum amplitude c at reference site) c associated temporal lag = (phase/360)*1/freq c c c determine reference phase -- in degrees c do j=1,nmode c c calculate magnitude and relative phase c do i=1,iabv angle(i)=avephase(i,j) cmag(i)=S(j)*aveamp(i,j)*envmax(j)*snorm(i)/weight(i) end do c c absolute phase at this point has no meaning c c need to rotate to reference phase convention c c complex pattern format c if (iref.eq.0) then angle_ref = 0.0 else angle_ref = angle(iref) endif if (iformat.eq.0) then do i=1,iabv angle(i)=angle(i)-angle_ref if (freq.gt.zero) then tlag = real(angle(i))/(360.0*freq) write (33+j-1,2345) i,real(cmag(i)),real(angle(i)), $ tlag else write (33+j-1,1234) i,real(cmag(i)),real(angle(i)) endif end do else c c real-valued successive snapshot format c phase_max = 180.0*(icycle+1) phase_step = phase_max/float(nsnap) do ii=1,nsnap+1 do i=1,iabv phase = float(ii-1)*phase_step-angle(i)+angle(iref) phase_ref = float(ii-1)*phase_step tlag = zero if (freq.gt.zero) then tlag = phase_ref/(360.0*freq) endif cr = cmag(i)*cos(phase*pi/180.0) write (33+j-1,2345) i,phase_ref,cr,tlag end do end do end if c c end do 1234 format (i4,2f14.4) 2345 format (i4,3f14.4) c c calculate total and partial fractional variances for each c mode c do j=1,nwin sum = 0.d0 do k=j,nwin sum = sum + S(k)**2 end do partvar(j)=0.0 if (sum.gt.1d-6) partvar(j)=S(j)**2/sum end do sum = 0.d0 do i=1,nwin sum = sum + S(i)**2 end do c c write (9,*) 'Signal reconstruction information:' write (9,*) write (9,*) '(average) frequency (cyc/year) =',freq if (freqsum.gt.zero) $ write (9,*) '(average) Period = ',one/freq,' year' write (9,*) 'mode sing. val frac var part var' do i=1,nwin write (9,999) i,real(S(i)),real(S(i)**2/sum),real(partvar(i)) end do 999 format (i2,f12.2,2f12.4) if (nsites.gt.0) then write (9,*) if (ilow.gt.0) then write (9,*) '% lowpassed variance explained' else write (9,*) '% raw variance explained' endif write (9,*) 'for mode # ',imoderecon,' at gridpoints:' sumtot1 = 0.0 sumtot2 = 0.0 do i=1,nsites iii = igrid(i) sum1(iii)=0.0 sum2(iii)=0.0 do im=n1,n2 rcon(iii,im)=rconbest(iii,im) term1 = ((rcon(iii,im)-alow(im,iii)) $ *weight(iii)/snorm(iii))**2 term2 = (alow(im,iii)*weight(iii)/snorm(iii))**2 sum1(iii) = sum1(iii)+term1 sumtot1 = sumtot1+term1 sum2(iii) = sum2(iii)+term2 sumtot2 = sumtot2+term2 end do write (9,987) iii, 100*(one-sum1(iii)/sum2(iii)) end do if (ilow.gt.0) then write (9,*) 'total % lowpassed variance over all ', $ nsites,' sites: ',100*(one-sumtot1/sumtot2) else write (9,*) 'total % raw variance over all ', $ nsites,' sites: ',100*(one-sumtot1/sumtot2) endif if (ilambda.eq.1) then write (6,*) 'error minimization:' write (6,*) ' lambda1 lambda2 lambda3 ', $ lambmin1,lambmin2,lambmin3 write (9,*) 'error minimization:' write (9,*) ' lambda1 lambda2 lambda3 ', $ lambmin1,lambmin2,lambmin3 if (ilow.gt.0) then write (6,*) ' % lowpassed var', $ 100.0*(one-summintot1/summintot2) write (9,*) ' % lowpassed var', $ 100.0*(one-summintot1/summintot2) else write (6,*) ' % raw var', $ 100.0*(one-summintot1/summintot2) write (9,*) ' % raw var', $ 100.0*(one-summintot1/summintot2) endif endif endif 986 format (i3,2f10.4,2x,f7.4,2x,f6.1) 987 format (i4,f6.1) c c close output files c close (unit=9) do i=1,nmode close (unit=13+i-1) close (unit=23+i-1) close (unit=33+i-1) end do write (6,*) 'done.' stop end c SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V) C C Singular value decomposition of an M by N complex matrix A, C where M .GT. N . The singular values are stored in the vector C S. The first NU columns of the M by M unitary matrix U and the C first NV columns of the N by N unitary matrix V that minimize C Det(A-USV*) are also computed. C C C P.A. Businger and G.H. Golub, "Singular Value Decomposition C of a Complex Matrix," Communications of the ACM, vol. 12, C pp. 564-565, October 1969. C C This algorithm is reprinted by permission, Association for C Computing Machinery; copyright 1969. C parameter (nwork=10000) COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R REAL *8 S(NMAX),B(nwork),C(nwork),T(nwork),zero,one,ETA,TOL, $ EPS ETA = 1.2d-7 TOL = 2.4d-32 zero = 0.d0 one = 1.d0 NP=N+IP N1=N+1 C Householder reduction C(1)=zero K=1 10 K1=K+1 C Elimination of A(I,K) , I=K+1,...,M Z=zero DO 20 I=K,M 20 Z=Z+REAL(A(I,K))**2+dIMAG(A(I,K))**2 B(K)=zero IF (Z .LE. TOL) GO TO 70 Z=SQRT(Z) B(K)=Z W=CdABS(A(K,K)) Q=dcmplx(one,zero) IF (W .NE. zero) Q=A(K,K)/W A(K,K)=Q*(Z+W) IF (K .EQ. NP) GO TO 70 DO 50 J=K1,NP Q=dcmplx(zero,zero) DO 30 I=K,M 30 Q=Q+dCONJG(A(I,K))*A(I,J) Q=Q/(Z*(Z+W)) DO 40 I=K,M 40 A(I,J)=A(I,J)-Q*A(I,K) 50 CONTINUE C Phase transformation Q=-dCONJG(A(K,K))/CdABS(A(K,K)) DO 60 J=K1,NP 60 A(K,J)=Q*A(K,J) C Elimination of A(K,J) , J=K+2,...,N 70 IF (K .EQ. N) GO TO 140 Z=zero DO 80 J=K1,N 80 Z=Z+REAL(A(K,J))**2+dIMAG(A(K,J))**2 C(K1)=zero IF (Z .LE. TOL) GO TO 130 Z=SQRT(Z) C(K1)=Z W=CdABS(A(K,K1)) Q=dcmplx(one,zero) IF (W .NE. zero) Q=A(K,K1)/W A(K,K1)=Q*(Z+W) DO 110 I=K1,M Q=dcmplx(zero,zero) DO 90 J=K1,N 90 Q=Q+dCONJG(A(K,J))*A(I,J) Q=Q/(Z*(Z+W)) DO 100 J=K1,N 100 A(I,J)=A(I,J)-Q*A(K,J) 110 CONTINUE C Phase transformation Q=-dCONJG(A(K,K1))/CdABS(A(K,K1)) DO 120 I=K1,M 120 A(I,K1)=A(I,K1)*Q 130 K=K1 GO TO 10 C Tolerance for negligible elements 140 EPS=zero DO 150 K=1,N S(K)=B(K) T(K)=C(K) 150 EPS=DMAX1(EPS,S(K)+T(K)) EPS=EPS*ETA C Initialization of U and V IF (NU .EQ. 0) GO TO 180 DO 170 J=1,NU DO 160 I=1,M 160 U(I,J)=dcmplx(zero,zero) 170 U(J,J)=dcmplx(one,zero) 180 IF (NV .EQ. 0) GO TO 210 DO 200 J=1,NV DO 190 I=1,N 190 V(I,J)=dcmplx(zero,zero) 200 V(J,J)=dcmplx(one,zero) C QR diagonalization 210 DO 380 KK=1,N K=N1-KK C Test for split 220 DO 230 LL=1,K L=K+1-LL IF (ABS(T(L)) .LE. EPS) GO TO 290 IF (ABS(S(L-1)) .LE. EPS) GO TO 240 230 CONTINUE C Cancellation of B(L) 240 CS=zero SN=one L1=L-1 DO 280 I=L,K F=SN*T(I) T(I)=CS*T(I) IF (ABS(F) .LE. EPS) GO TO 290 H=S(I) W=SQRT(F*F+H*H) S(I)=W CS=H/W SN=-F/W IF (NU .EQ. 0) GO TO 260 DO 250 J=1,N X=REAL(U(J,L1)) Y=REAL(U(J,I)) U(J,L1)=dCMPLX(X*CS+Y*SN,0.) 250 U(J,I)=dCMPLX(Y*CS-X*SN,0.) 260 IF (NP .EQ. N) GO TO 280 DO 270 J=N1,NP Q=A(L1,J) R=A(I,J) A(L1,J)=Q*CS+R*SN 270 A(I,J)=R*CS-Q*SN 280 CONTINUE C Test for convergence 290 W=S(K) IF (L .EQ. K) GO TO 360 C Origin shift X=S(L) Y=S(K-1) G=T(K-1) H=T(K) F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y) G=SQRT(F*F+one) IF (F .LT. zero) G=-G F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X C QR step CS=one SN=one L1=L+1 DO 350 I=L1,K G=T(I) Y=S(I) H=SN*G G=CS*G W=SQRT(H*H+F*F) T(I-1)=W CS=F/W SN=H/W F=X*CS+G*SN G=G*CS-X*SN H=Y*SN Y=Y*CS IF (NV .EQ. 0) GO TO 310 DO 300 J=1,N X=REAL(V(J,I-1)) W=REAL(V(J,I)) V(J,I-1)=dCMPLX(X*CS+W*SN,0.) 300 V(J,I)=dCMPLX(W*CS-X*SN,0.) 310 W=SQRT(H*H+F*F) S(I-1)=W CS=F/W SN=H/W F=CS*G+SN*Y X=CS*Y-SN*G IF (NU .EQ. 0) GO TO 330 DO 320 J=1,N Y=REAL(U(J,I-1)) W=REAL(U(J,I)) U(J,I-1)=dCMPLX(Y*CS+W*SN,0.) 320 U(J,I)=dCMPLX(W*CS-Y*SN,0.) 330 IF (N .EQ. NP) GO TO 350 DO 340 J=N1,NP Q=A(I-1,J) R=A(I,J) A(I-1,J)=Q*CS+R*SN 340 A(I,J)=R*CS-Q*SN 350 CONTINUE T(L)=zero T(K)=F S(K)=X GO TO 220 C Convergence 360 IF (W .GE. zero) GO TO 380 S(K)=-W IF (NV .EQ. 0) GO TO 380 DO 370 J=1,N 370 V(J,K)=-V(J,K) 380 CONTINUE C Sort singular values DO 450 K=1,N G=-one J=K DO 390 I=K,N IF (S(I) .LE. G) GO TO 390 G=S(I) J=I 390 CONTINUE IF (J .EQ. K) GO TO 450 S(J)=S(K) S(K)=G IF (NV .EQ. 0) GO TO 410 DO 400 I=1,N Q=V(I,J) V(I,J)=V(I,K) 400 V(I,K)=Q 410 IF (NU .EQ. 0) GO TO 430 DO 420 I=1,N Q=U(I,J) U(I,J)=U(I,K) 420 U(I,K)=Q 430 IF (N .EQ. NP) GO TO 450 DO 440 I=N1,NP Q=A(J,I) A(J,I)=A(K,I) 440 A(K,I)=Q 450 CONTINUE C Back transformation IF (NU .EQ. 0) GO TO 510 DO 500 KK=1,N K=N1-KK IF (B(K) .EQ. zero) GO TO 500 Q=-A(K,K)/CdABS(A(K,K)) DO 460 J=1,NU 460 U(K,J)=Q*U(K,J) DO 490 J=1,NU Q=dcmplx(zero,zero) DO 470 I=K,M 470 Q=Q+dCONJG(A(I,K))*U(I,J) Q=Q/(CdABS(A(K,K))*B(K)) DO 480 I=K,M 480 U(I,J)=U(I,J)-Q*A(I,K) 490 CONTINUE 500 CONTINUE 510 IF (NV .EQ. 0) GO TO 570 IF (N .LT. 2) GO TO 570 DO 560 KK=2,N K=N1-KK K1=K+1 IF (C(K1) .EQ. zero) GO TO 560 Q=-dCONJG(A(K,K1))/CdABS(A(K,K1)) DO 520 J=1,NV 520 V(K1,J)=Q*V(K1,J) DO 550 J=1,NV Q=dcmplx(zero,zero) DO 530 I=K1,N 530 Q=Q+A(K,I)*V(I,J) Q=Q/(CdABS(A(K,K1))*C(K1)) DO 540 I=K1,N 540 V(I,J)=V(I,J)-Q*dCONJG(A(K,I)) 550 CONTINUE 560 CONTINUE 570 RETURN END c c subroutine lowpass(xx,yy,nscan,fcut,dt) c c c filters from Stearns and David (c) c c Michael Mann c parameter (zero=0.0,one=1.0,two=2.0,M=4096) integer j,nscan real xx(nscan),yy(nscan) real x(0:M),g(M+1),b(0:M),c(0:M),f1,f2 c11 write (6,*) 'type of filter' c write (6,*) '(1) low-pass' c write (6,*) '(2) high-pass' c write (6,*) '(3) bandpass' c write (6,*) '(4) bandstop' c read (5,*) iband iband = 1 c22 write (6,*) 'type of filter window' c write (6,*) '(1) rectangular' c write (6,*) '(2) tapered rectangular' c write (6,*) '(3) triangular' c write (6,*) '(4) Hanning' c write (6,*) '(5) Hamming' c write (6,*) '(6) Blackman' c read (5,*) iwndo iwndo = 5 c c sampling periods in years c tau = dt c c normalized frequency is c tau * frequency(yr-1) c c write (6,*) 'frequencies in (yr^-1)' c write (6,*) 'cutoff frequency' f1 = fcut f2 = 0.0 f1 = f1*tau f2 = f2*tau c do i2=0,nscan-1 g(i2+1)=float(i2+1) x(i2)=xx(i2+1) b(i2)=zero end do N = i2 LL = N-1 call spfird(N-1,iband,f1,f2,iwndo,b,ierror) if (ierror.ne.0) then write (6,*) 'invalid frequencies' stop endif do i=0,2*N c(i)=zero end do do i=0,N-1 do j=0,N-1 c(i+j)=c(i+j)+x(i)*b(j) end do end do do i=0,N-1 yy(i+1)=c(i+N/2) end do return end c- subroutine spconv(x,y,l,nmin,ierror) c-latest date: 11/27/85 c-fast convolution of sequences x and y. c-x and y should be different vectors, dimensioned x(0:l) and y(0:l). c-l=last index in both x and y. must be (power of 2)+1 and at least 5. c-nmin=minimum shift of interest in the convolution function. c-fft length ,n, used internally, is l-1. c-let k=index of last nonzero sample in y(0)---y(n-1). then x(0) c- thru x(n-1) must include padding of at least k-nmin-1 zeros. c-convolution function (output) replaces x(nmin) thru x(n-1). c-ierror=0 no error detected c- 1 l-1 not a power of 2 c- 2 nmin out of range c- 3 inadequate zero padding dimension x(0:l),y(0:l) complex cx ierror=1 n=l-1 test=n 1 test=test/2. if(test-2.) 8,2,1 2 ierror=2 if(nmin.lt.0.or.nmin.ge.n) return ierror=3 do 3 k=n-1,0,-1 if(y(k).ne.0.) go to 4 3 continue 4 do 5 j=n-1,0,-1 if(x(j).ne.0.) go to 6 5 continue 6 if(n-j.le.k-nmin) return call spfftr(x,n) call spfftr(y,n) do 7 m=0,n/2 cx=cmplx(x(2*m),x(2*m+1))*cmplx(y(2*m),y(2*m+1)) x(2*m)=real(cx)/n x(2*m+1)=aimag(cx)/n 7 continue call spiftr(x,n) ierror=0 8 return end c- subroutine spfftr(x,n) c-latest date: 11/12/85 c-fft routine for real time series (x) with n=2**k samples. c-computation is in place, output replaces input. c-dimension x(0:n+1) (real; not complex) or larger. c-real time series (input) is in x(0) thru x(n-1). c-first output fft component is x(0)=real part; x(1)=imaginary, c-etc. last component is x(n)=real and x(n+1)=imaginary part. c-important: n must be at least 4 and must be a power of 2. complex x(0:n/2),u,tmp tpn=8.*atan(1.)/n call spfftc(x,n/2,-1) x(n/2)=x(0) do 1 m=0,n/4 u=cmplx(sin(m*tpn),cos(m*tpn)) tmp=((1.+u)*x(m)+(1.-u)*conjg(x(n/2-m)))/2. x(m)=((1.-u)*x(m)+(1.+u)*conjg(x(n/2-m)))/2. x(n/2-m)=conjg(tmp) 1 continue return end c- subroutine spiftr(x,n) c-latest date: 02/20/87 c-inverse fft of the complex spectrum of a real time series. c-x and n are the same as in spfftr. important: n must be a power c-of 2 and x must be dimensioned x(0:n+1) (real array, not complex). c-this routine transforms the output of spfftr back into the input, c-scaled by n. computation is in place, as in spfftr. complex x(0:n/2),u1,tmp tpn=8.*atan(1.)/n do 1 m=0,n/4 u1=cmplx(sin(m*tpn),-cos(m*tpn)) tmp=(1.+u1)*x(m)+(1.-u1)*conjg(x(n/2-m)) x(m)=(1.-u1)*x(m)+(1.+u1)*conjg(x(n/2-m)) x(n/2-m)=conjg(tmp) 1 continue call spfftc(x,n/2,1) return end c- subroutine spfftc(x,n,isign) c-latest date: 02/20/87 c-fast fourier transform of n=2**k complex data points using time c-decomposition with input bit reversal. n must be a power of 2. c-x must be specified complex x(0:n-1)or larger. c-input is n complex samples, x(0),x(1),...,x(n-1). c-computation is in place, output replaces input. c-isign = -1 for forward transform, +1 for inverse. c-x(0) becomes the zero transform component, x(1) the first, c-and so forth. x(n-1) becomes the last component. complex x(0:n-1),t pisign=4*isign*atan(1.) mr=0 do 2 m=1,n-1 l=n 1 l=l/2 if(mr+l.ge.n) go to 1 mr=mod(mr,l)+l if(mr.le.m) go to 2 t=x(m) x(m)=x(mr) x(mr)=t 2 continue l=1 3 if(l.ge.n) return do 5 m=0,l-1 do 4 i=m,n-1,2*l t=x(i+l)*exp(cmplx(0.,m*pisign/float(l))) x(i+l)=x(i)-t x(i)=x(i)+t 4 continue 5 continue l=2*l go to 3 end c function spwndo(itype,n,k) c-latest date: 11/13/85 c-this function generates a single sample of a data window. c-itype=1(rectangular), 2(tapered rectangular), 3(triangular), c- 4(hanning), 5(hamming), or 6(blackman). c- (note: tapered rectangular has cosine-tapered 10% ends.) c-n=size (total no. samples) of window. c-k=sample number within window, from 0 through n-1. c- (if k is outside this range, spwndo is set to 0.) pi=4.*atan(1.) spwndo=0. if(itype.lt.1.or.itype.gt.6) return if(k.lt.0.or.k.ge.n) return spwndo=1. go to (1,2,3,4,5,6), itype 1 return 2 l=(n-2)/10 if(k.le.l) spwndo=0.5*(1.0-cos(k*pi/(l+1))) if(k.gt.n-l-2) spwndo=0.5*(1.0-cos((n-k-1)*pi/(l+1))) return 3 spwndo=1.0-abs(1.0-2*k/(n-1.0)) return 4 spwndo=0.5*(1.0-cos(2*k*pi/(n-1))) return 5 spwndo=0.54-0.46*cos(2*k*pi/(n-1)) return 6 spwndo=0.42-0.5*cos(2*k*pi/(n-1))+0.08*cos(4*k*pi/(n-1)) return end c c subroutine spfird(l,iband,fln,fhn,iwndo,b,ierror) c-latest date: 05/18/86 c-fir digital filter design using windowed fourier series c-l=length of filter = l+1 c-iband=1(lowpass); 2(highpass); 3(bandpass); 4(bandstop) c-fln=normalized low cut-off frequency in hz-sec c-fhn=normalized high cut-off (bp,bs) in hz-sec c-iwndo=1(rectangular); 2(tapered rectangular); 3(triangular) c- 4(hanning); 5(hamming); 6(blackman) c-digital filter coefficients are returned in b(0:l) c-ierror=0 no errors detected c- 1 invalid length (l<=0) c- 2 invalid window type c- 3 invalid filter type c- 4 invalid cut-off frequency dimension b(0:l) pi=4.*atan(1.) do 1 i=0,l b(i)=0. 1 continue ierror=0 if(l.le.0) ierror=1 if(iwndo.lt.1.or.iwndo.gt.6) ierror=2 if(iband.lt.1.or.iband.gt.4) ierror=3 if(fln.le.0..or.fln.gt.0.5) ierror=4 if(iband.ge.3.and.fln.ge.fhn) ierror=4 if(iband.ge.3.and.fhn.ge.0.5) ierror=4 if(ierror.ne.0) return dly=l/2. lim=int(l/2) mid=0 if(dly.eq.lim) then lim=lim-1 mid=1 endif wc1=2.*pi*fln if(iband.ge.3) wc2=2.*pi*fhn go to (5,10,15,20) iband c-lowpass design 5 continue do 6 i=0,lim s=i-dly b(i)=((sin(wc1*s))/(pi*s))*spwndo(iwndo,l+1,i) b(l-i)=b(i) 6 continue if(mid.eq.1) b(l/2)=wc1/pi return c-highpass design 10 continue do 11 i=0,lim s=i-dly b(i)=((sin(pi*s)-sin(wc1*s))/(pi*s))*spwndo(iwndo,l+1,i) b(l-i)=b(i) 11 continue if(mid.eq.1) b(l/2)=1.-wc1/pi return c-bandpass design 15 continue do 16 i=0,lim s=i-dly b(i)=((sin(wc2*s)-sin(wc1*s))/(pi*s))*spwndo(iwndo,l+1,i) b(l-i)=b(i) 16 continue if(mid.eq.1) b(l/2)=(wc2-wc1)/pi return c-bandstop design 20 continue do 21 i=0,lim s=i-dly b(i)=((sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s))* + spwndo(iwndo,l+1,i) b(l-i)=b(i) 21 continue if(mid.eq.1) b(l/2)=(wc1+pi-wc2)/pi return end c c subroutine taper(n,nwin,el) c c to generate slepian tapers c ta is a real*4 array c c j. park c real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb real*8 dfac,drat,gamma,bh,ell common/nnaamm/iwflag common/npiprol/anpi common/tapsum/tapsum(20),ell(20) common/work/ip(8194) common/taperzz/z(65536) ! we use this common block for scratch space common/stap2/ta(8192,16) dimension a(8192,8),el(10) data iwflag/0/,pi/3.14159265358979d0/ equivalence (a(1,1),ta(1,1)) an=dfloat(n) ww=dble(anpi)/an cs=dcos(2.d0*pi*ww) c initialize matrix for eispack subroutine c type *,'ww,cs,an',ww,cs,an do i=0,n-1 ai=dfloat(i) a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2 a(i+1,2)=-ai*(an-ai)/2.d0 a(i+1,3)=a(i+1,2)**2 ! required by eispack routine end do eps=1.e-13 m11=1 call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip, x ierr,a(1,4),a(1,5)) c type *,ierr,rlb,rlu c call tnou('eigenvalues for the eigentapers') c type *,(el(i),i=1,nwin) call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr, x a(1,4),a(1,5),a(1,6),a(1,7),a(1,8)) c we calculate the eigenvalues of the dirichlet-kernel problem c i.e. the bandwidth retention factors c from slepian 1978 asymptotic formula, gotten from thomson 1982 eq 2.5 c supplemented by the asymptotic formula for k near 2n from slepian 1978 eq 61 dfac=an*pi*ww drat=8.d0*dfac dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac) do k=1,nwin el(k)=1.d0-dfac dfac=dfac*drat/k ! is this correct formula? yes,but fails as k -> 2n end do c call tnou('eigenvalues for the eigentapers (small k)') c type *,(el(i),i=1,nwin) gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0 do k=1,nwin bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma ell(k)=1.d0/(1.d0+dexp(pi*bh)) end do c call tnou('eigenvalues for the eigentapers (k near 2n)') c type *,(ell(i),i=1,nwin) do i=1,nwin el(i)=dmax1(ell(i),el(i)) end do c call tnou('composite asymptotics for taper eigenvalues') c type *,(el(i),i=1,nwin) c normalize the eigentapers to preserve power for a white process c i.e. they have rms value unity do k=1,nwin kk=(k-1)*n tapsum(k)=0. tapsq=0. do i=1,n aa=z(kk+i) ta(i,k)=aa tapsum(k)=tapsum(k)+aa tapsq=tapsq+aa*aa end do aa=sqrt(tapsq/n) tapsum(k)=tapsum(k)/aa do i=1,n ta(i,k)=ta(i,k)/aa end do end do c type *,'tapsum',(tapsum(i),i=1,nwin) c refft preserves amplitudes with zeropadding c zum beispiel: a(i)=1.,i=1,100 will transform at zero freq to b(f=0.)=100 c no matter how much zero padding is done c therefore we need not doctor the taper normalization, c but wait until the fft to force the desired units iwflag=1 return end