subroutine FFTEXE c C Intermediate routine which prepares input for the FFT calculation c on the recorded interferogram. C C The amount of zero-filling is determined by the variable MULZER - this c can be set externally to this routine (but given a negative sign), in which c case no screen output is made by this routine. C PARAMETER (Nmaxpt=8192,maxpts=400) c real data(nmaxpt),w1(2*nmaxpt),w2(nmaxpt) real p(nmaxpt) integer idata(maxpts) real*8 fstep logical ovrlap common /points/data,npts common /work/w1 common /work1/w2 common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /specf/p,fstep,ncut,MULZER,NCALL C c c...Transfer points to array DATA passed to FFT routines: c minimum number of points is 256 or higher powers of 2 so fill-up c appropriately c NCALL=NCALL+1 NOWRIT=0 i=0 do 170 j=nskips+1,nrep-nskipe i=i+1 data(i)=idata(j) 170 continue npts=i IF(NPTS.LE.8192)NMAX=8192 if(npts.le.4096)nmax=4096 if(npts.le.2048)nmax=2048 if(npts.le.1024)nmax=1024 if(npts.le.512)nmax=512 if(npts.le.256)nmax=256 if(npts.ne.nmax)then do 52 i=npts+1,nmax data(i)=data(npts) 52 continue npts=nmax endif c c...FFT options: amount of zero-filling: note that -ve MULZER can be C generated by routine LOOKIN for recalculation of FFT with previous c value of MULZER, ie. when cutoffs for rejected points have been modified C C -on first time into this routine (NCALL=1) default value of MULZER is c assigned c IF(MULZER.GE.0)THEN write(*,56)npts 56 format(1x/i7,' points readied for FFT') 31 if(NCALL.ne.1)then 33 write(*,30) 30 format(1x/' Specify zero filling parameter n defined as ' * ' Total_Npts = 2**n * Npts'/ * ' (n=0,1,2,3.... )'// * 45x,' n = ',$) read(*,*,err=31)mulzer if(npts*2**mulzer.lt.npts)goto 33 else mulzer=4 555 if(npts*2**mulzer.gt.nmaxpt)then MULZER=MULZER-1 GOTO 555 ENDIF ENDIF C if(npts*2**mulzer.gt.nmaxpt)then mulzer=npts*2**mulzer write(*,775)mulzer,nmaxpt,char(7) 775 format(1x/' ERROR:',i8,' points needed for FFT but only ', * i6,' dimensioned',a// * ' -----> Specify smaller n'/) goto 31 endif ELSE mulzer=-mulzer 556 if(npts*2**mulzer.gt.nmaxpt)then MULZER=MULZER-1 GOTO 556 ENDIF NOWRIT=1 ENDIF C SHIFT=DATA(NPTS) do 32 i=1,npts*2**mulzer IF(I.LE.NPTS)THEN DATA(I)=DATA(I)-SHIFT ELSE DATA(I)=0.0 ENDIF 32 continue npts=npts*2**mulzer C k=1 m=npts/(k+1) IF(NOWRIT.NE.1)write(*,57)npts 57 format(1x/i7,' point record will be used for FFT') c c...FFT c samfre=1./tstep fnyq=samfre/2. fstep=samfre/npts ncut=fnyq/fstep fstep=fstep/1000. IF(NOWRIT.NE.1) * write(*,25)tstep*1.E+6,samfre*1.E-6,fnyq*1.E-6,fstep,ncut 25 format(1x/' time step = ',f15.10,' microsec.'/ * ' sampling frequency = ',f15.10,' MHz'/ * ' Nyquist frequency = ',f15.10,' MHz'/ * ' frequency step = ',f15.10,' kHz'/ * ' points in Nyq. interval = ',i15/) c ovrlap=.true. call spctrm(m,k,ovrlap) C C...Postscale points in power spectrum so that their intensities for c various amounts of zero filling are unified and equivalent to those c for n=4 c ymult=mulzer-4 ymult=16.d0**nint(ymult) if(nint(ymult).ne.1)then do 570 i=1,ncut p(i)=p(i)*ymult 570 continue endif c return end C C------------------------------------------------------------------------ c SUBROUTINE spctrm(m,k,ovrlap) c c Power spectrum estimation using routine 'four1' c c p = on output contains the input data's power (mean square amplitude) at c frequency (j-1)/(2*m) cycles per gridpoint, for j=1,2....,m c m = number of data points in segment c k = number of segments (each with 2m data points) c ovrlap=.false. segments do not overlap, 4*m*k data points c ovrlap=.true. segments overlap, (2k+1)*m data points c data = time domain data points c parameter (nmax=8192) INTEGER k,m c REAL p(m),w1(2*nmax),w2(nmax) REAL p(nmax),w1(2*nmax),w2(nmax) LOGICAL ovrlap real data(nmax) real*8 fstep common /points/data,npts common /work/w1 common /work1/w2 common /specf/p,fstep,ncut,NFFT,NCALL c INTEGER j,j2,joff,joffn,kk,m4,m43,m44,mm REAL den,facm,facp,sumw,w,window window(j)=(1.-abs(((j-1)-facm)*facp)) Bartlett c window(j)=1. Square c window(j)=(1.-(((j-1)-facm)*facp)**2 Welch nread=0 mm=m+m m4=mm+mm m44=m4+4 m43=m4+3 den=0. facm=m facp=1./m c c...accumulate the squared sum of the weights c sumw=0. do 11 j=1,mm sumw=sumw+window(j)**2 11 continue c c...initialize the spectrum to zero c do 12 j=1,m p(j)=0. 12 continue c c...initialize the 'save' half-buffer - this is a modifcation to use c the data from common block /points/. The values are read in c successively and NREAD is the total number of data points used. c If more points are required then are in the data then the last point c is repeated c if(ovrlap)then do 21 j=1,m nread=nread+1 if(nread.gt.npts)then w2(j)=data(npts) else w2(j)=data(nread) endif 21 continue endif c c...Loop over data set segments in groups of two. Get two complete c segments into workspace. c do 18 kk=1,k do 15 joff=-1,0,1 if(ovrlap)then do 13 j=1,m w1(joff+j+j)=w2(j) 13 continue do 22 j=1,m nread=nread+1 if(nread.gt.npts)then w2(j)=0. else w2(j)=data(nread) endif 22 continue joffn=joff+mm do 14 j=1,m w1(joffn+j+j)=w2(j) 14 continue else do 23 j=joff+2,m4,2 nread=nread+1 if(nread.gt.npts)then w1(j)=0. else w1(j)=data(nread) endif 23 continue endif 15 continue c c...Apply the window to the data c do 16 j=1,mm j2=j+j w=window(j) w1(j2)=w1(j2)*w w1(j2-1)=w1(j2-1)*w 16 continue c c...Fourier transform the windowed data c call four1(mm,1) c c...Sum results into previous segments c p(1)=p(1)+w1(1)**2+w1(2)**2 do 17 j=2,m j2=j+j p(j)=p(j)+w1(j2)**2+w1(j2-1)**2 * +w1(m44-j2)**2+w1(m43-j2)**2 17 continue den=den+sumw 18 continue c c...Correct normalization and normalize the output c den=m4*den do 19 j=1,m p(j)=p(j)/den 19 continue c write(*,25)nread 25 format(1x/i10,' points used in FFT') c return end c c---------------------------------------------------------------------------- c SUBROUTINE four1(nn,isign) parameter (nmax=8192) INTEGER isign,nn c REAL data(2*nn) real data(2*nmax) common /work/data c c Routine replaces data(1:2*nn) by its discrete Fourier transform, if isign c is input as -1; or replaces data(1:2*nn) by nn times its inverse discrete c Fourier transform, if isign is input as -1. c data is a complex array of length nn, or equivalently, a real array of c length 2*nn. c nn MUST be an integer power of 2 (this is not checked for!) c INTEGER i,istep,j,m,mmax,n REAL tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp n=2*nn j=1 do 11 i=1,n,2 if(j.gt.i)then tempr=data(j) tempi=data(j+1) data(j)=data(i) data(j+1)=data(i+1) data(i)=tempr data(i+1)=tempi endif m=n/2 1 if((m.ge.2).and.(j.gt.m))then j=j-m m=m/2 goto 1 endif j=j+m 11 continue c mmax=2 2 if(n.gt.mmax)then istep=2*mmax theta=6.28318530717959d0/(isign*mmax) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 do 13 m=1,mmax,2 do 12 i=m,n,istep j=i+mmax tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1) tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j) data(j)=data(i)-tempr data(j+1)=data(i+1)-tempi data(i)=data(i)+tempr data(i+1)=data(i+1)+tempi 12 continue wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi 13 continue mmax=istep goto 2 endif c return end c C------------------------------------------------------------------------ c SUBROUTINE PF(fmark,smax,XPEAK,ERRORX,htcut,xbot,fhalm,fhalp) C C Position of line maximum is established using the bisection method: C midpoints of line contour are determined at preselected number of sections c NSECT from peak maximum down to a selected fraction of peak height, HTCUT. c C Straight line fit to such midpoints gives, for y equal to line maximum, c the line frequency, with some account for possible line asymmetry. C C FMARK - frequency of the marker, which is assumed to have been set near C peak maximum c SMAX - on exit the intensity at the maximum C XPEAK - on exit the required central fitted peak position (requires C addition of FCENT) C ERRORX - on exit error on the fitted peak position c HTCUT - proportion of line height to which profile division is taken c XBOT - value of X for Y=0.4Ymax for drawing fitted line c FHALM - offset for negative FWHH point c FHALP - offset for positive FWHH point C PARAMETER (Nmaxpt=8192,nsect=20) c real p(nmaxpt),fr(nsect,2) real*8 fstep,fmark REAL*8 SUMX,SUMY,SUMXY,SUMX2,SUMY2,CXX,CXY,CYY,RN,A0,A1, * xoffs,yoffs c common /specf/p,fstep,npts,NFFT,NCALL common /peak/x(nsect),y(nsect) C c...determine initial value of peak maximum and its position c n=(fmark/fstep)+1. SMAX=P(N) 1 IF(N.LE.1.OR.N.GE.NPTS-1)THEN XPEAK=(N-1)*FSTEP ERRORX=(NPTS-1)*FSTEP RETURN ENDIF IF(P(N+1).GT.SMAX)THEN SMAX=P(N+1) N=N+1 GOTO 1 ENDIF IF(P(N-1).GT.SMAX)THEN SMAX=P(N-1) N=N-1 GOTO 1 ENDIF NMAX=N C C...Determine frequencies of points on sections through line contour: c linear interpolation used C do 2 ns=1,nsect ysect=smax-ns*smax*htcut/real(nsect) y(ns)=ysect do 3 n=nmax,npts if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fr(ns,2)=(rn-1.0)*fstep goto 2 endif 3 continue 2 continue c do 5 ns=1,nsect ysect=y(ns) do 4 n=nmax,2,-1 if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fr(ns,1)=(rn-1.0)*fstep goto 5 endif 4 continue 5 continue C do 6 n=1,nsect x(n)=0.5*(fr(n,1)+fr(n,2)) 6 continue c C...Straight line fit: since the line is almost vertical (gradient very c large) this was found to lead to numerical instabilities and for c this reason axes are reversed for least squares. For further increased c numerical stability SMAX is subtracted from Y and X(1) from X ie. the c equation of fit is: c c (x-x1) = a0 + a1 (y-smax) C SUMX=0.D0 SUMY=0.D0 SUMXY=0.D0 SUMX2=0.D0 SUMY2=0.D0 XOFFS=X(1) YOFFS=SMAX DO 7 I=1,NSECT SUMy=SUMy+(X(I)-xoffs) SUMx=SUMx+(Y(I)-yoffs) SUMXY=SUMXY+(X(I)-xoffs)*(Y(I)-yoffs) SUMy2=SUMy2+(X(I)-xoffs)**2 SUMx2=SUMx2+(Y(I)-yoffs)**2 7 CONTINUE C C...coefficients RN=NSECT CXX=SUMX2-SUMX*SUMX/RN CXY=SUMXY-SUMX*SUMY/RN A1=CXY/CXX IF(A1.EQ.0.D0)THEN ERRORX=0.D0 XPEAK=0.D0 RETURN ENDIF A0=(SUMY-A1*SUMX)/RN c c...peak frequency xpeak=a0+x(1) c C...coordinates of 0.4Ymax point xbot=a0-0.6*smax*a1+x(1) C C...error CYY=SUMY2-SUMY*SUMY/RN ERA1S=((CYY/CXX)-(CXY/CXX)**2)/(RN-2.D0) ERA0S=SUMX2*ERA1S/RN ERRORX=dsqrt(dble(era0s)) C C...Find X values at FWHH C ysect=0.5*smax do 13 n=nmax,npts if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fhalp=(rn-1.0)*fstep goto 12 endif 13 continue fhalp=(npts-1)*fstep c 12 do 14 n=nmax,2,-1 if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fhalm=(rn-1.0)*fstep goto 15 endif 14 continue fhalm=fstep 15 continue c RETURN END C C------------------------------------------------------------------------ C SUBROUTINE SORTH c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c parameter (maxspe=70,maxpts=400,nivols=7) character fnams(maxspe)*12,ftemp(maxspe)*12,filarc*30 real detvol(maxspe,2),volint(maxspe,nivols) integer interf(maxspe,maxpts) integer*2 ipt(maxspe),iseen(maxspe) real*8 wk(maxspe) common /scans/interf,wk,ipt,filarc common /scans2/detvol,volint,fnams,iseen,nscans c INTEGER*2 IIPT,L,N,NSTART,I,J,IR REAL*8 WWK real rtemp(maxspe) equivalence (rtemp(1),ftemp(1)) C nstart=1 n=iabs(nscans) c L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT GOTO 100 ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c c...reorder FNAMS, DETVOL and VOLINT according to the order in IPT c (this is not done if NSCANS has previously been made negative) c 100 if(nscans.lt.0)return do 101 l=1,2 do 102 i=1,nscans rtemp(i)=detvol(i,l) 102 continue do 101 i=1,nscans j=ipt(i) detvol(i,l)=rtemp(j) 101 continue c do 201 l=1,nivols do 202 i=1,nscans rtemp(i)=volint(i,l) 202 continue do 201 i=1,nscans j=ipt(i) volint(i,l)=rtemp(j) 201 continue c do 302 i=1,nscans ftemp(i)=fnams(i) 302 continue do 301 i=1,nscans j=ipt(i) fnams(i)=ftemp(j) 301 continue c c RETURN END C C_____________________________________________________________________________ c subroutine baksub(nspt) c PARAMETER (Nmaxpt=8192,maxpts=400,maxsmo=199) c integer idata(maxpts),ioldat(maxpts),itemp(maxpts) real spol(maxsmo) common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /smooth/ioldat,itemp,spol c C...Subtraction of background by triple smoothing of interferogram with c least squares smoothing interval c C For smoothing interval of length 2m+1 the elements of the smoothing c (cubic) polynomial are given by: C C 3(3m**2 + 3m -1 - 5s**2) C c(s) = ------------------------ C (2m+1) (2m-1) (2m+3) C C where s runs from -m to +m (T.H.Edwards and P.D.Wilson, Applied C Spectroscopy 28,541-545(1974)) C C c...set up coefficients in smoothing polynomial M=NSPT/2 T1=3.D0/((2*M+1)*(2.D0*M-1.D0)*(2*M+3)) T2=3*M*M+3.D0*M-1.D0 DO 1103 j=1,NSPT IS=j-M-1 SPOL(j)=T1*(T2-5*IS*IS) 1103 CONTINUE C c...Smooth three times ISTRT=M+1 IFIN=nrep-M DO 1104 k=1,3 do 1105 j=1,nrep itemp(j)=idata(j) 1105 continue DO 543 I=1,nrep SUM=0. DO 544 J=1,NSPT IS=J-M-1 II=I+IS IF(II.LT.1)II=-II+1 IF(II.GT.nrep)II=nrep-(II-nrep-1) SUM=SUM+itemp(II)*SPOL(J) 544 CONTINUE Idata(I)=sum 543 CONTINUE DO 545 I=1,nrep ITEMP(I)=IDATA(I) 545 CONTINUE 1104 continue C do 1106 j=1,nrep idata(j)=ioldat(j)-idata(j) 1106 continue c return end C C_____________________________________________________________________________ C_____________________________________________________________________________