Ticket #934: smoo.f

File smoo.f, 17.1 KB (added by Roman Tolchenov, 11 years ago)
Line 
1C SRSLIB PROGRAM SMFFT
2C                     
3C PURPOSE: SMOOTHS DATA BY FILTERING OUT HIGH HARMONICS IN FT SPECTRUM
4C             
5C LIMITATIONS:
6C             
7C 1. ARRAYS DIMENSIONED FOR MAXIMUM 32800 POINTS
8C 2. NO CHECKS FOR VALID INPUT DATA ARE MADE   
9C 3. ABSCISSA VALUES ASSUMED EQUIDISTANT       
10C 4. NAG FFT ROUTINES 
11C
12C GLOBALS:
13C
14C CDSQ   : LOGARITHM OF POWER SPECTRUM 
15C WF     : WIENER FILTER ORDINATE (DECIBELS)
16C RNOISE : NOISE LEVEL                     
17C RJLEV  : REJECTION LEVEL                 
18C ITURN  : NUMBER OF ITERATIONS             
19C NUMBER : NO OF OUTLIER POINTS             
20C J0     : THESHOLD FREQUENCY               
21C
22C LOCAL VARIABLES: 
23C                 
24C X     : ABSCISSA VALUES
25C YORIG : ORIGINAL ORDINATE VALUES 
26C YANS  : SMOOTHED ORDINATE VALUES 
27C YDIF  : RESIDUAL ORDINATE VALUES 
28C TITLE : INPUT FILE TITLE RECORD 
29C ANSWER: USER ANSWER TO PROMPT   
30C IDSN  : INPUT FILE NAME         
31C IREJ  : OUTLIER REJECTION FLAG
32C FOUR  : CONSTANT 4.3429448   
33C EPS7  : CONSTANT 10E-7       
34C IPOS  : SCREEN POSITION POINTER
35C NPOINT: NUMBER OF INPUT POINTS 
36C                               
37        SUBROUTINE SMOO(smoo_get, smoo_put)
38        EXTERNAL smoo_get, smoo_put
39        PARAMETER (np=512,nph=256,mn=33000)
40        COMMON /WIEN/ CDSQ(mn),WF(mn),RNOISE,ITURN,NUMBER,J0,RJLEV   
41        real*4 Xin(mn),Yin(mn),Ein(mn),Xout(mn),Yout(mn),Eout(mn)
42        LOGICAL*1 IDSN(44)     
43        real*4 X(mn), YORIG(mn), YANS(mn), YDIF(mn)     
44
45        DATA IREJ/1/,IPOS/770/
46        DATA FOUR/4.3429448/,EPS7/10.0E-7/
47        DATA INPUT/1/,ITERM/5/,IPRINT/6/,IDEV/5/,IOUT/17/,IWF/18/       
48C
49C=========READ INPUT DATA                                               00000970
50C                                                                       00000980
51        call module_get_int(smoo_get, 'lptin', lptin)   
52        call module_get_real_array(smoo_get, 'Xin', Xin, lptin) 
53        call module_get_real_array(smoo_get, 'Yin', Yin, lptin) 
54        call module_get_real_array(smoo_get, 'Ein', Ein, lptin) 
55        call module_get_int(smoo_get, 'ider', ider)     
56
57        if(lptin.eq.0)then              !check #points
58          call module_error(" smoo>",
59        1 "ERROR  ** No input data", " ")
60          call module_put_int(smoo_put, 'kill', 1)     
61          RETURN
62        endif
63        do n=1,lptin
64         X(n)=xin(n)
65         yorig(n)=yin(n)
66        end do
67
68        npoint= lptin
69        irej = 1
70C         
71C=========FILTER DATA
72C                   
73        jflag=0
74c       ider=0
75      CALL FILTER (X,YORIG,YANS,YDIF,NPOINT,IREJ,NHALF,JFLAG,IDER)           
76
77        do I=1,NHALF
78         IF (CDSQ(I).LT.EPS7) CDSQ(I)=EPS7
79         IF (WF(I).LT.EPS7) WF(I)=EPS7   
80         CDSQ(I)=FOUR*ALOG(CDSQ(I)/RNOISE)
81         WF(I)=FOUR*ALOG(WF(I))           
82        end do             
83
84        lptout=lptin
85        do n=1,lptout
86         yout(n)=yans(n)
87        end do
88        call module_put_int(smoo_put, 'lptout', LPTOUT)       
89        call module_put_real_array(smoo_put, 'Xout', Xout, lptout)     
90        call module_put_real_array(smoo_put, 'Yout', Yout, lptout)     
91        call module_put_real_array(smoo_put, 'Eout', Eout, lptout)
92        call module_put_int(smoo_put, 'kill', 0)       
93c
94        RETURN
95        END
96C++++++++++++++++++++++++++++++++++++++++++++++++
97      SUBROUTINE FILTER (X,YOR,Y,YDIF,NPOINT,IREJ,NHALF,JFLAG,IDER)   
98C                                                     
99C PURPOSE: FOURIER FILTERS INPUT DATA                 
100C                                                     
101        parameter (mn=33000)
102        real*4 X(mn),YOR(mn),Y(mn),YDIF(mn),Z(mn)
103        real*4 DER1(mn),DER2(mn)
104        INTEGER NPOINT,IREJ,NHALF,JFLAG
105C                                                       
106C INPUT ARGUMENTS:                                       
107C                                                       
108C     X       : ABSCISSA VALUES                         
109C     YOR     : ORDINATE VALUES                           
110C     NPOINT  : NUMBER OF COORDINATE VALUES             
111C     JFLAG   : CUT-OFF FREQUENCY OVERRIDING J0 IF JFLAG>0   
112C     IREJ    : OUTLIER REJECTION FLAG, 2 FOR YES, 1 FOR NO   
113C                                                 
114C RETURNED VALUES :                               
115C                                               
116C     Y       : SMOOTHED ORDINATES             
117C     YDIF    : RESIDUALS                       
118C     NHALF   : NUMBER OF FOURIER SPECTRUM POINTS   
119C                                       
120C GLOBALS:                               
121C                                             
122        COMMON /WIEN/ CDSQ(mn),WF(mn),RNOISE,ITURN,NUMBER,J0,RJLEV 
123        REAL*4 CDSQ,WF,RNOISE,RJLEV
124        INTEGER ITURN,NUMBER,J0
125C                                                   
126C CALLS   5: REJECT , C06EAF , FILTER , C06GBF , C06EBF 
127C CALLED BY: MAIN                                 
128C                                                 
129C LOCAL VARIABLES:                               
130C                                               
131        DOUBLE PRECISION COEFF(mn),A(mn),B(mn),C(mn),D(mn)
132        dimension LIER(mn),ILST(26)
133        INTEGER N2,M,MPLUS1,M1,NFULL,N0,KK,K,N0K1
134        LOGICAL INVERS
135C                                                 
136C A, B  :REAL AND AND IMAGINARY PART OF FT     
137C LIERS : INDICES OF OUTLIER POINTS           
138C                                             
139        DATA ZERO/0.0/,P6/0.6931451/,IPRINT/6/
140C                                               
141        ITURN=0                                     
142        ICOUNT = IDER
143        NUMBER=0
144C                                       
145C=========FIND OUTLIER REJECTION LEVEL FOR NPOINT POINTS
146C                                               
147C      CALL REJECT (NPOINT,RJLEV)               
148C      WRITE (IPRINT,250) NPOINT,RJLEV           
149C                                             
150C=========FIND M POWER OF 2 SUCH THAT 2**M<N<=2**(M+1)   
151C                                             
152        IFLAG = 0
153        M=ALOG(FLOAT(NPOINT))/P6
154        IF (NPOINT.EQ.2**M) M=M-1
155        MPLUS1=M+1
156        NHALF=2**M
157C                                           
158C=========FIND DIFFERENCE BETWEEN NPOINT AND 2**MPLUS1 
159C                                   
160        NFULL=2*NHALF
161        MDIFF=NFULL-NPOINT
162        WRITE (IPRINT,280) MDIFF,MPLUS1,NFULL
163        do I=1,NPOINT
164         Y(I)=YOR(I)
165        end do
166C                                 
167C=========SUBTRACT CUBIC SPLINE TIED AT FIRST AND LAST POINTS 
168C                                         
169 20     continue
170        ICOUNT = ICOUNT - 1
171
172        CALL SUBSPLINE(NPOINT,X,Y,Z,DER1,DER2)
173      DO 30 I=1,NPOINT
174         YI=Y(I)
175         Y(I)=YI-Z(I)
176  30  CONTINUE
177C                                 
178C=========SET POINTS BETWEEN NPOINT AND 2**MPLUS1 TO 0 IF MDIFF>0 
179C                                     
180      IF (MDIFF.EQ.0) GOTO 50
181      DO 40 I=1,MDIFF
182         Y(NPOINT+I)=ZERO
183  40  CONTINUE
184C                                           
185C=========FOURIER TRANSFORM                           
186C         VALUES RETURNED IN A AND B ARE THE           
187C         THE COSINE AND SINE COEFFICIENTS, RESPECTIVELY. 
188C                                               
189  50  IFAIL=0
190        ITURN=ITURN+1
191        IF ( IFLAG .EQ. 0 ) THEN
192         do I=1,NFULL
193          COEFF(I)=Y(I)
194         end do
195        ELSE
196         JFLAG = 0
197         do I=1,NFULL
198          COEFF(I) = Z(I)
199         end do
200        ENDIF
201
202        CALL C06EAF (COEFF,NFULL,IFAIL)
203c      WRITE (IPRINT,200) NFULL
204C                                     
205C=========ASSIGN FOURIER COEFFICIENTS     
206C                                           
207      N2=(NFULL+1)/2
208c      WRITE (IPRINT,210) N2
209      A(1)=COEFF(1)
210      B(1)=0.0
211      DO 70 J=2,N2
212         NJ=NFULL-J+2
213         A(J)=COEFF(J)
214         A(NJ)=COEFF(J)
215         B(J)=COEFF(NJ)
216         B(NJ)=-COEFF(NJ)
217  70  CONTINUE
218      IF (MOD(NFULL,2).NE.0) GOTO 80
219      A(N2+1)=COEFF(N2+1)
220      B(N2+1)=0.0
221C                                           
222C=========CALCULATE POWER SPECTRUM CDSQ   
223C                                 
224      FPOINT=FLOAT(NFULL)/4.0
225      N21=N2+1
226  80  DO 90 I=1,N21
227         AA=A(I)
228         BB=B(I)
229         CDSQ(I)=(AA*AA+BB*BB)/FPOINT
230  90  CONTINUE
231      FFFF=FLOAT(NFULL)/FLOAT(NPOINT)
232      NPHALF=NPOINT/2
233      K=NPHALF/2
234      SUMN0=ZERO
235      FNK=FLOAT(NPHALF)/(FLOAT(2*K))
236      DO 100 KK=1,K
237         J=FLOAT(N21)-FLOAT(KK)*FFFF+1.0
238         SUMN0=SUMN0+CDSQ(J)
239  100 CONTINUE
240C                                               
241C=========FIND FLUCTUATION LEVEL OF FOURIER HARMONICS 
242C                                           
243        RNOISE=SUMN0/FLOAT(2*K)
244        SIGMA2=FNK*SUMN0
245        SIGMA=SQRT(SIGMA2)
246c      WRITE (IPRINT,290) SIGMA,RNOISE
247C                                       
248C=========CALCULATE FILTER FUNCTION AND MULTIPLY FOURIER COEFFICIENTS
249C                                       
250        CALL WIENER (A,B,N21,JFLAG)
251c      WRITE (IPRINT,220) N2,J0       
252        write(IPRINT,220) N2
253C                                     
254C=========INVERSE FOURIER TRANSFORM   
255C                                 
256        IFLAG = IFLAG + 1
257        ISTEM= IDER/2
258        STEM= -1*ISTEM
259        IF (IDER.EQ.0 .OR. ICOUNT.LT.IDER - 1) THEN
260
261       DO 110 J=2,N2
262         C(J)=A(J)
263         D(J)=B(J)
264         NJ=NFULL-J+2
265         C(NJ)=C(J)
266         D(NJ)=-D(J)
267
268  110  CONTINUE
269
270      ELSE IF (IDER.EQ.1) THEN
271
272       DO 111 J=2,N2
273         C(J)=FLOAT(J)*-B(J)
274         D(J)=FLOAT(J)*A(J)
275         NJ=NFULL-J+2
276         C(NJ)=C(J)
277         D(NJ)=-D(J)
278  111 CONTINUE
279
280      ELSE IF ( IDER .EQ. 2 ) THEN
281
282       DO 112 J=2,N2
283         CTEM= STEM*FLOAT(J)**IDER
284         C(J)=CTEM*A(J)
285         D(J)=CTEM*B(J)
286         NJ=NFULL-J+2
287         C(NJ)=C(J)
288         D(NJ)=-D(J)
289  112  CONTINUE
290      ENDIF
291
292
293      COEFF(1)=A(1)
294      DO 120 J=2,N2
295         NJ=NFULL-J+2
296         COEFF(J)=0.5*(C(J)+C(NJ))
297         COEFF(NJ)=0.5*(D(J)-D(NJ))
298  120 CONTINUE
299      IF (MOD(NFULL,2).EQ.0) COEFF(N2+1)=A(N2+1)
300C                                         
301        CALL C06GBF (COEFF,NFULL,IFAIL)
302        CALL C06EBF (COEFF,NFULL,IFAIL)
303        IDERLAST= IDER
304c      WRITE (IPRINT,230) NFULL                 
305C                                             
306C=========RESTORE ORDINATES TO ORIGINAL SCALE BY ADDING "BACKGROUND"
307C                                           
308      IF (IFLAG .EQ. 1 ) THEN
309
310         DO I=1,NPOINT
311           Y(I) = COEFF(I)
312         END DO
313         IF ( IDER .EQ. 0 ) THEN
314         DO I=1,NPOINT
315            Y(I) = Y(I) + Z(I)
316         END DO
317         GOTO 9999
318         ENDIF
319      ELSE
320         DO I=1,NPOINT
321           Z(I) = COEFF(I)
322           Y(I) = Y(I) + Z(I)
323         END DO
324      ENDIF
325
326      IF ( IFLAG .LT. 2 ) GOTO 50
327C
328C      IF (IDER .EQ. 0 .OR. ICOUNT .LT. IDER - 1 ) THEN
329C       DO 130 I=1,NPOINT             
330C         YI=COEFF(I)               
331C         TYPE *,YI
332C         Y(I)=YI+Z(I)
333C  130  CONTINUE                     
334C      ELSE IF (IDER .EQ. 1) THEN
335C       DO 131 I=1,NPOINT           
336C         YI=COEFF(I)             
337C         Y(I)=YI+DER1(I)
338C  131  CONTINUE                     
339C      ELSE IF (IDER .EQ. 2) THEN
340C       DO 132 I=1,NPOINT             
341C         YI=COEFF(I)               
342C         Y(I)=YI+DER2(I)
343C  132  CONTINUE                     
344C      END IF
345c9999   WRITE (IPRINT,240) NPOINT   
3469999    continue
347C      JFLAG = 0
348C      J0 = 0
349C      IF ( ICOUNT .GT. -1 ) GOTO 20
350C                                   
351C=========CALCULATE RESIDUALS     
352C                                   
353      DO 140 I=1,NPOINT                 
354         YDIF(I)=YOR(I)-Y(I)
355  140 CONTINUE
356      IF (IREJ.EQ.1) GOTO 180
357C                                           
358C=========FIND OUTLIERS FOR CYCLE ITURN     
359C                                           
360      IDN=NUMBER
361      NUMBER=0
362      DO 150 I=1,NPOINT
363         DIF=ABS(YDIF(I))/SIGMA
364         IF (DIF.LE.RJLEV) GOTO 150
365         NUMBER=NUMBER+1
366         LIER(NUMBER)=I
367  150 CONTINUE
368c      WRITE (IPRINT,300) NUMBER,ITURN       
369c      WRITE (IPRINT,260) (LIER(K),K=1,NUMBER) 
370        IF (NUMBER.LE.0) GOTO 180
371        IDN=NUMBER-IDN
372        IF (IDN.EQ.0) GOTO 180
373C                                     
374C=========REPLACE OUTLIER ORDINATE WITH CURRENT SMOOTHED ORDINATE VALUE
375C                                                                     
376c      WRITE (IPRINT,310)             
377      DO 170 I=1,NPOINT
378         DO 160 K=1,NUMBER
379            IF (I.EQ.LIER(K)) GOTO 170
380  160    CONTINUE
381         Y(I)=YOR(I)
382  170 CONTINUE
383      GOTO 20
384C                                           
385C=========NO MORE OUTLIERS TO BE REPLACED 
386C                                       
387c  180 WRITE (IPRINT,270) NUMBER,ITURN   
388180     continue
389        RETURN
390C                             
391  200 FORMAT(' > FORWARD FFT FOR ',I5,' POINTS COMPLETE')
392  210 FORMAT(' > ASSIGNING ',I5,' FOURIER COEFFICIENTS')
393c  220 FORMAT(' > WIENER FILTER CALCULATED FOR '
394c     + ,I5,' COEFFICIENTS',/,
395c     + 11x,' CUT-OFF FREQUENCY ',I4)
396220     FORMAT(' fourier> filter at ',I5,' coefficients')
397  230 FORMAT(' > INVERSE FFT FOR ',I5,' POINTS COMPLETE')
398  240 FORMAT(' > BASELINE ORDINATES ADDED TO ',I5,' POINTS')
399  250 FORMAT (11x,' REJECTION LEVEL FOR ',I5,' POINTS:',E13.5)
400  260 FORMAT (128(8I5,/))
401  270 FORMAT (11x,I6,' OUTLIER POINTS REJECTED IN ',I3,' CYCLES')
402  280 FORMAT (' fourier> ',I6,' points set to zero to padd 2**',I4,
403        1'=',I5,' data points')
404  290 FORMAT (11x,' SIGMA=',E13.5,' RNOISE=',E13.5)
405  300 FORMAT (11x,I6,' OUTLIER POINTS FOR CYCLE',I5)
406  310 FORMAT (' REPEATING FILTERING TO REJECT OUTLIERS')
407        END
408
409C+++++++++++++++++++++++++++++++
410        SUBROUTINE REJECT (N,X)
411C                                             
412C PURPOSE: COMPUTATES OUTLIER REJECTION LEVEL FOR RANDOM NUMBERS 
413C          WITH NORMAL DISTRIBUTION       
414C                           
415C CALLS   0:                 
416C CALLED BY: FILTER         
417C                           
418        NP=2*(N/2)
419        X=3.8+0.15*(ALOG(FLOAT(N))/ALOG(2.0)-6.0)
420C                                                           
421C...USE NEXT LINE OR IBM SCIENTIFIC SUBROUTINE LIBRARY ROUTINE NDTRI 
422C                                       
423CIBM  P=(1.0-0.01/2.0)**(1.0/FLOAT(NP))
424CIBM      CALL NDTRI (P,X,D,IER) 
425        RETURN
426        END
427
428C++++++++++++++++
429        SUBROUTINE WIENER (A,B,N,JFLAG)
430C                               
431C PURPOSE: CALCULATES WIENER FILTER FUNCTION 
432C                                 
433        parameter (mn=33000)
434        DOUBLE PRECISION A(mn),B(mn)
435        INTEGER N,JFLAG
436C                               
437        COMMON /WIEN/ CDSQ(mn),WF(mn),RNOISE,ITURN,NUMBER,J0,RJLEV
438        DOUBLE PRECISION C,D,WFI,ALPHA
439        REAL*4 CDSQ,WF,RNOISE,RJLEV
440        INTEGER ITURN,NUMBER,J0
441C                                                       
442C     WF= CD1/(1+CD1),  WHERE CD1=CD/RNOISE, AND  CD=C*C+D*D   
443C     AS  CD1>1                                         
444C     AND EXTRAPOLATION OF CD1 UP TO J1, WHERE J1 EQUALS       
445C     THE HARMONIC NUMBER, WHERE CD1 (DB) < -20 DB.     
446C                                                       
447C CALLS   0:                                           
448C CALLED BY: FILTER                                     
449C                                                       
450              DATA ONE/1.0/,TWO/2.0/,ZERO/0.0/,FOUR/4.3429448/,DBLEV/-20.0/
451C                                               
452c      WRITE(6,1)
453c1     FORMAT(/,11X,'>> INPUT ALPHA :'$)
454c      ACCEPT *,ALPHA
455        do IK=1,N
456         WF(IK)=ZERO
457        end do
458        J0=1
459        J1=1
460        XX=ZERO
461        XY=ZERO
462        YM=ZERO
463      DO 40 J=1,N
464         CD1=CDSQ(J)/RNOISE
465         CD2=ALOG(CD1)
466         IF (JFLAG.NE.0) GOTO 20
467         IF (J0.GT.1) GOTO 50
468         IF ((CD1.LE.ONE).AND.(J0.EQ.1)) J0=J
469         GOTO 30
470  20     IF (JFLAG.LT.J) GOTO 50
471  30    WF(J)=CD1/(ONE+CD1)
472         XX=XX+FLOAT(J*J)
473         XY=XY+FLOAT(J)*CD2
474         YM=YM+CD2
475  40  CONTINUE
476  50  IF (JFLAG.NE.0) J0=JFLAG
477        RJ0F=FLOAT(J0)
478        XM=(ONE+RJ0F)/TWO
479        YM=YM/RJ0F
480        A1=(XY-RJ0F*XM*YM)/(XX-RJ0F*XM**2)
481        B1=YM-A1*XM
482        J1=(DBLEV/FOUR-B1)/A1
483        do J=J0,J1
484         S=EXP(A1*FLOAT(J)+B1)
485         WF(J)=S/(ONE+S)
486        end do
487c     ALPHA = 0.2
488c     
489        do I=1,N
490         WFI=WF(I)
491         C=A(I)
492         D=B(I)
493c        WFI = WFI / ( 1. + DEXP ( DFLOAT ( I - J0 ) / ( ALPHA *
494c    +       DFLOAT ( J0 ) ) ) )
495         A(I)=WFI*C
496         B(I)=WFI*D
497        end do
498        RETURN
499        END
500
501        subroutine SUBSPLINE(n,x,y,z,der1,der2)
502c
503c..     This subroutine evaluates a crude background z(i)
504c       to be used as a baseline subtraction
505c       and also 1st to 3rd derivatives.
506        parameter  (nmax=33000,knots=10)
507        real*4  x(n),y(n),z(nmax),der1(nmax),der2(nmax)
508        real*4  val(4)
509        real*4  xknot(knots),yknot(knots)
510        real*4  dknot(knots),work(3*knots)
511c
512        ntem= n/(knots-1)
513        xknot(1)= x(1)
514        xknot(knots)= x(n)
515        yknot(1)= y(1)
516        yknot(knots)= y(n)
517c
518        do l=1,knots-2
519c       iimin= l*ntem - jnint(0.5*float(ntem))
520         iimin= l*ntem - nint(0.5*float(ntem))
521         call MINMAX2(y(iimin),ntem,yknot(l+1),yt,imin,imax)
522         item= iimin + imin -1
523         xknot(l+1)= x(item)
524        end do
525c
526        call TB04A(knots,xknot,yknot,dknot,work)
527c
528        iopt= 1
529        do i= 1,n
530         call TG02A(iopt,knots,xknot,yknot,dknot,x(i),val)
531         z(i)= val(1)
532         der1(i)= val(2)
533         der2(i)= val(3)
534        end do
535c
536        return
537        end
538c
539        subroutine MINMAX2(z,n,zmin,zmax,imin,imax)
540        real*4 z(n)
541c
542c This subroutine finds the minimum and maximum values in a graph
543c and additionally gives the bin no. associated with these values.
544c
545        zmin=z(1)
546        zmax=z(1)
547        imin=1
548        imax=1
549        do 10 i=1,n
550         if (z(i) .gt. zmax) then
551                zmax=z(i)
552                imax=i
553                goto 10
554         end if
555         if (z(i) .lt. zmin) then
556                zmin=z(i)
557                imin=i
558         end if
559 10     continue
560c
561        return
562        end
563* CMS REPLACEMENT HISTORY, Element MINMAX2.FOR
564* *1    13-APR-1988 10:45:13 CMM "V2.3 Baseline elements inserted"
565* CMS REPLACEMENT HISTORY, Element MINMAX2.FOR