1 | C SRSLIB PROGRAM SMFFT |
---|
2 | C |
---|
3 | C PURPOSE: SMOOTHS DATA BY FILTERING OUT HIGH HARMONICS IN FT SPECTRUM |
---|
4 | C |
---|
5 | C LIMITATIONS: |
---|
6 | C |
---|
7 | C 1. ARRAYS DIMENSIONED FOR MAXIMUM 32800 POINTS |
---|
8 | C 2. NO CHECKS FOR VALID INPUT DATA ARE MADE |
---|
9 | C 3. ABSCISSA VALUES ASSUMED EQUIDISTANT |
---|
10 | C 4. NAG FFT ROUTINES |
---|
11 | C |
---|
12 | C GLOBALS: |
---|
13 | C |
---|
14 | C CDSQ : LOGARITHM OF POWER SPECTRUM |
---|
15 | C WF : WIENER FILTER ORDINATE (DECIBELS) |
---|
16 | C RNOISE : NOISE LEVEL |
---|
17 | C RJLEV : REJECTION LEVEL |
---|
18 | C ITURN : NUMBER OF ITERATIONS |
---|
19 | C NUMBER : NO OF OUTLIER POINTS |
---|
20 | C J0 : THESHOLD FREQUENCY |
---|
21 | C |
---|
22 | C LOCAL VARIABLES: |
---|
23 | C |
---|
24 | C X : ABSCISSA VALUES |
---|
25 | C YORIG : ORIGINAL ORDINATE VALUES |
---|
26 | C YANS : SMOOTHED ORDINATE VALUES |
---|
27 | C YDIF : RESIDUAL ORDINATE VALUES |
---|
28 | C TITLE : INPUT FILE TITLE RECORD |
---|
29 | C ANSWER: USER ANSWER TO PROMPT |
---|
30 | C IDSN : INPUT FILE NAME |
---|
31 | C IREJ : OUTLIER REJECTION FLAG |
---|
32 | C FOUR : CONSTANT 4.3429448 |
---|
33 | C EPS7 : CONSTANT 10E-7 |
---|
34 | C IPOS : SCREEN POSITION POINTER |
---|
35 | C NPOINT: NUMBER OF INPUT POINTS |
---|
36 | C |
---|
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/ |
---|
48 | C |
---|
49 | C=========READ INPUT DATA 00000970 |
---|
50 | C 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 |
---|
70 | C |
---|
71 | C=========FILTER DATA |
---|
72 | C |
---|
73 | jflag=0 |
---|
74 | c 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) |
---|
93 | c |
---|
94 | RETURN |
---|
95 | END |
---|
96 | C++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
97 | SUBROUTINE FILTER (X,YOR,Y,YDIF,NPOINT,IREJ,NHALF,JFLAG,IDER) |
---|
98 | C |
---|
99 | C PURPOSE: FOURIER FILTERS INPUT DATA |
---|
100 | C |
---|
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 |
---|
105 | C |
---|
106 | C INPUT ARGUMENTS: |
---|
107 | C |
---|
108 | C X : ABSCISSA VALUES |
---|
109 | C YOR : ORDINATE VALUES |
---|
110 | C NPOINT : NUMBER OF COORDINATE VALUES |
---|
111 | C JFLAG : CUT-OFF FREQUENCY OVERRIDING J0 IF JFLAG>0 |
---|
112 | C IREJ : OUTLIER REJECTION FLAG, 2 FOR YES, 1 FOR NO |
---|
113 | C |
---|
114 | C RETURNED VALUES : |
---|
115 | C |
---|
116 | C Y : SMOOTHED ORDINATES |
---|
117 | C YDIF : RESIDUALS |
---|
118 | C NHALF : NUMBER OF FOURIER SPECTRUM POINTS |
---|
119 | C |
---|
120 | C GLOBALS: |
---|
121 | C |
---|
122 | COMMON /WIEN/ CDSQ(mn),WF(mn),RNOISE,ITURN,NUMBER,J0,RJLEV |
---|
123 | REAL*4 CDSQ,WF,RNOISE,RJLEV |
---|
124 | INTEGER ITURN,NUMBER,J0 |
---|
125 | C |
---|
126 | C CALLS 5: REJECT , C06EAF , FILTER , C06GBF , C06EBF |
---|
127 | C CALLED BY: MAIN |
---|
128 | C |
---|
129 | C LOCAL VARIABLES: |
---|
130 | C |
---|
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 |
---|
135 | C |
---|
136 | C A, B :REAL AND AND IMAGINARY PART OF FT |
---|
137 | C LIERS : INDICES OF OUTLIER POINTS |
---|
138 | C |
---|
139 | DATA ZERO/0.0/,P6/0.6931451/,IPRINT/6/ |
---|
140 | C |
---|
141 | ITURN=0 |
---|
142 | ICOUNT = IDER |
---|
143 | NUMBER=0 |
---|
144 | C |
---|
145 | C=========FIND OUTLIER REJECTION LEVEL FOR NPOINT POINTS |
---|
146 | C |
---|
147 | C CALL REJECT (NPOINT,RJLEV) |
---|
148 | C WRITE (IPRINT,250) NPOINT,RJLEV |
---|
149 | C |
---|
150 | C=========FIND M POWER OF 2 SUCH THAT 2**M<N<=2**(M+1) |
---|
151 | C |
---|
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 |
---|
157 | C |
---|
158 | C=========FIND DIFFERENCE BETWEEN NPOINT AND 2**MPLUS1 |
---|
159 | C |
---|
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 |
---|
166 | C |
---|
167 | C=========SUBTRACT CUBIC SPLINE TIED AT FIRST AND LAST POINTS |
---|
168 | C |
---|
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 |
---|
177 | C |
---|
178 | C=========SET POINTS BETWEEN NPOINT AND 2**MPLUS1 TO 0 IF MDIFF>0 |
---|
179 | C |
---|
180 | IF (MDIFF.EQ.0) GOTO 50 |
---|
181 | DO 40 I=1,MDIFF |
---|
182 | Y(NPOINT+I)=ZERO |
---|
183 | 40 CONTINUE |
---|
184 | C |
---|
185 | C=========FOURIER TRANSFORM |
---|
186 | C VALUES RETURNED IN A AND B ARE THE |
---|
187 | C THE COSINE AND SINE COEFFICIENTS, RESPECTIVELY. |
---|
188 | C |
---|
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) |
---|
203 | c WRITE (IPRINT,200) NFULL |
---|
204 | C |
---|
205 | C=========ASSIGN FOURIER COEFFICIENTS |
---|
206 | C |
---|
207 | N2=(NFULL+1)/2 |
---|
208 | c 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 |
---|
221 | C |
---|
222 | C=========CALCULATE POWER SPECTRUM CDSQ |
---|
223 | C |
---|
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 |
---|
240 | C |
---|
241 | C=========FIND FLUCTUATION LEVEL OF FOURIER HARMONICS |
---|
242 | C |
---|
243 | RNOISE=SUMN0/FLOAT(2*K) |
---|
244 | SIGMA2=FNK*SUMN0 |
---|
245 | SIGMA=SQRT(SIGMA2) |
---|
246 | c WRITE (IPRINT,290) SIGMA,RNOISE |
---|
247 | C |
---|
248 | C=========CALCULATE FILTER FUNCTION AND MULTIPLY FOURIER COEFFICIENTS |
---|
249 | C |
---|
250 | CALL WIENER (A,B,N21,JFLAG) |
---|
251 | c WRITE (IPRINT,220) N2,J0 |
---|
252 | write(IPRINT,220) N2 |
---|
253 | C |
---|
254 | C=========INVERSE FOURIER TRANSFORM |
---|
255 | C |
---|
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) |
---|
300 | C |
---|
301 | CALL C06GBF (COEFF,NFULL,IFAIL) |
---|
302 | CALL C06EBF (COEFF,NFULL,IFAIL) |
---|
303 | IDERLAST= IDER |
---|
304 | c WRITE (IPRINT,230) NFULL |
---|
305 | C |
---|
306 | C=========RESTORE ORDINATES TO ORIGINAL SCALE BY ADDING "BACKGROUND" |
---|
307 | C |
---|
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 |
---|
327 | C |
---|
328 | C IF (IDER .EQ. 0 .OR. ICOUNT .LT. IDER - 1 ) THEN |
---|
329 | C DO 130 I=1,NPOINT |
---|
330 | C YI=COEFF(I) |
---|
331 | C TYPE *,YI |
---|
332 | C Y(I)=YI+Z(I) |
---|
333 | C 130 CONTINUE |
---|
334 | C ELSE IF (IDER .EQ. 1) THEN |
---|
335 | C DO 131 I=1,NPOINT |
---|
336 | C YI=COEFF(I) |
---|
337 | C Y(I)=YI+DER1(I) |
---|
338 | C 131 CONTINUE |
---|
339 | C ELSE IF (IDER .EQ. 2) THEN |
---|
340 | C DO 132 I=1,NPOINT |
---|
341 | C YI=COEFF(I) |
---|
342 | C Y(I)=YI+DER2(I) |
---|
343 | C 132 CONTINUE |
---|
344 | C END IF |
---|
345 | c9999 WRITE (IPRINT,240) NPOINT |
---|
346 | 9999 continue |
---|
347 | C JFLAG = 0 |
---|
348 | C J0 = 0 |
---|
349 | C IF ( ICOUNT .GT. -1 ) GOTO 20 |
---|
350 | C |
---|
351 | C=========CALCULATE RESIDUALS |
---|
352 | C |
---|
353 | DO 140 I=1,NPOINT |
---|
354 | YDIF(I)=YOR(I)-Y(I) |
---|
355 | 140 CONTINUE |
---|
356 | IF (IREJ.EQ.1) GOTO 180 |
---|
357 | C |
---|
358 | C=========FIND OUTLIERS FOR CYCLE ITURN |
---|
359 | C |
---|
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 |
---|
368 | c WRITE (IPRINT,300) NUMBER,ITURN |
---|
369 | c 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 |
---|
373 | C |
---|
374 | C=========REPLACE OUTLIER ORDINATE WITH CURRENT SMOOTHED ORDINATE VALUE |
---|
375 | C |
---|
376 | c 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 |
---|
384 | C |
---|
385 | C=========NO MORE OUTLIERS TO BE REPLACED |
---|
386 | C |
---|
387 | c 180 WRITE (IPRINT,270) NUMBER,ITURN |
---|
388 | 180 continue |
---|
389 | RETURN |
---|
390 | C |
---|
391 | 200 FORMAT(' > FORWARD FFT FOR ',I5,' POINTS COMPLETE') |
---|
392 | 210 FORMAT(' > ASSIGNING ',I5,' FOURIER COEFFICIENTS') |
---|
393 | c 220 FORMAT(' > WIENER FILTER CALCULATED FOR ' |
---|
394 | c + ,I5,' COEFFICIENTS',/, |
---|
395 | c + 11x,' CUT-OFF FREQUENCY ',I4) |
---|
396 | 220 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 | |
---|
409 | C+++++++++++++++++++++++++++++++ |
---|
410 | SUBROUTINE REJECT (N,X) |
---|
411 | C |
---|
412 | C PURPOSE: COMPUTATES OUTLIER REJECTION LEVEL FOR RANDOM NUMBERS |
---|
413 | C WITH NORMAL DISTRIBUTION |
---|
414 | C |
---|
415 | C CALLS 0: |
---|
416 | C CALLED BY: FILTER |
---|
417 | C |
---|
418 | NP=2*(N/2) |
---|
419 | X=3.8+0.15*(ALOG(FLOAT(N))/ALOG(2.0)-6.0) |
---|
420 | C |
---|
421 | C...USE NEXT LINE OR IBM SCIENTIFIC SUBROUTINE LIBRARY ROUTINE NDTRI |
---|
422 | C |
---|
423 | CIBM P=(1.0-0.01/2.0)**(1.0/FLOAT(NP)) |
---|
424 | CIBM CALL NDTRI (P,X,D,IER) |
---|
425 | RETURN |
---|
426 | END |
---|
427 | |
---|
428 | C++++++++++++++++ |
---|
429 | SUBROUTINE WIENER (A,B,N,JFLAG) |
---|
430 | C |
---|
431 | C PURPOSE: CALCULATES WIENER FILTER FUNCTION |
---|
432 | C |
---|
433 | parameter (mn=33000) |
---|
434 | DOUBLE PRECISION A(mn),B(mn) |
---|
435 | INTEGER N,JFLAG |
---|
436 | C |
---|
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 |
---|
441 | C |
---|
442 | C WF= CD1/(1+CD1), WHERE CD1=CD/RNOISE, AND CD=C*C+D*D |
---|
443 | C AS CD1>1 |
---|
444 | C AND EXTRAPOLATION OF CD1 UP TO J1, WHERE J1 EQUALS |
---|
445 | C THE HARMONIC NUMBER, WHERE CD1 (DB) < -20 DB. |
---|
446 | C |
---|
447 | C CALLS 0: |
---|
448 | C CALLED BY: FILTER |
---|
449 | C |
---|
450 | DATA ONE/1.0/,TWO/2.0/,ZERO/0.0/,FOUR/4.3429448/,DBLEV/-20.0/ |
---|
451 | C |
---|
452 | c WRITE(6,1) |
---|
453 | c1 FORMAT(/,11X,'>> INPUT ALPHA :'$) |
---|
454 | c 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 |
---|
487 | c ALPHA = 0.2 |
---|
488 | c |
---|
489 | do I=1,N |
---|
490 | WFI=WF(I) |
---|
491 | C=A(I) |
---|
492 | D=B(I) |
---|
493 | c WFI = WFI / ( 1. + DEXP ( DFLOAT ( I - J0 ) / ( ALPHA * |
---|
494 | c + 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) |
---|
502 | c |
---|
503 | c.. This subroutine evaluates a crude background z(i) |
---|
504 | c to be used as a baseline subtraction |
---|
505 | c 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) |
---|
511 | c |
---|
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) |
---|
517 | c |
---|
518 | do l=1,knots-2 |
---|
519 | c 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 |
---|
525 | c |
---|
526 | call TB04A(knots,xknot,yknot,dknot,work) |
---|
527 | c |
---|
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 |
---|
535 | c |
---|
536 | return |
---|
537 | end |
---|
538 | c |
---|
539 | subroutine MINMAX2(z,n,zmin,zmax,imin,imax) |
---|
540 | real*4 z(n) |
---|
541 | c |
---|
542 | c This subroutine finds the minimum and maximum values in a graph |
---|
543 | c and additionally gives the bin no. associated with these values. |
---|
544 | c |
---|
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 |
---|
560 | c |
---|
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 |
---|