'******************************************************************** '* SMOOTHING AN ARRAY OF N ORDINATES Y's (ASCENDING ORDER ABCISSAS) * '* ---------------------------------------------------------------- * '* Description: * '* This program uses the procedure SMOOFT for smoothing an array of * '* given ordinates (y's) that are in order of increasing abscissas * '* (x's), but without using the abscissas themselves supposed to be * '* equally spaced. It first removes any linear trend, then uses a * '* Fast Fourier Transform procedure (REALFT) to low-pass filter the * '* data. The linear trend is reinserted at the end. One user-speci- * '* fied parameter, EPS, enters "the amount of smoothing", given as * '* the number of points over which the data should be smoothed. * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* Input data file contains: * '* 1024 * '* 0.00000000000000E+0000 7.50000000000000E-0001 * '* 9.21288168207468E-0003 7.77368637910513E-0001 * '* 1.84257633641494E-0002 8.34466556277221E-0001 * '* 2.76386450462383E-0002 9.03071871110114E-0001 * '* 3.68515267282987E-0002 9.92958153417021E-0001 * '* 4.60644084103592E-0002 1.09195646826811E+0000 * '* 5.52772900924197E-0002 1.15230452277865E+0000 * '* 6.44901717745370E-0002 1.06763022290215E+0000 * '* 7.37030534565974E-0002 1.34541171127239E+0000 * '* 8.29159351386579E-0002 1.48611048393104E+0000 * '* 9.21288168207184E-0002 1.09349703210864E+0000 * '* 1.01341698502779E-0001 1.72386602840743E+0000 * '* 1.10554580184839E-0001 1.14317464708984E+0000 * '* ---------------------- ---------------------- * '* 9.37871355209791E+0000 2.43969819122867E+0001 * '* 9.38792643377383E+0000 2.42468007203424E+0001 * '* 9.39713931544975E+0000 2.42436619192304E+0001 * '* 9.40635219712567E+0000 2.42829449073179E+0001 * '* 9.41556507880159E+0000 2.42980085689633E+0001 * '* 9.42477796047751E+0000 2.43119449022633E+0001 * '* * '* Output file contains (here EPS=30): * '* * '* Time Y Smoothed Y * '* ---------------------------------------- * '* 0.000000 0.750000 0.785805 * '* 0.009213 0.777369 0.812665 * '* 0.018426 0.834467 0.840955 * '* 0.027639 0.903072 0.871028 * '* 0.036852 0.992958 0.903131 * '* 0.046064 1.091956 0.937348 * '* 0.055277 1.152305 0.973577 * '* 0.064490 1.067630 1.011530 * '* 0.073703 1.345412 1.050779 * '* 0.082916 1.486110 1.090820 * '* 0.092129 1.093497 1.131156 * '* 0.101342 1.723866 1.171384 * '* 0.110555 1.143175 1.211274 * '* -------- -------- -------- * '* 9.378709 24.396982 24.225096 * '* 9.387921 24.246801 24.247604 * '* 9.397134 24.243662 24.271116 * '* 9.406346 24.282946 24.295259 * '* 9.415559 24.298008 24.319866 * '* 9.424771 24.311945 24.344969 * '* * '* ---------------------------------------------------------------- * '* Reference: "Numerical Recipes By W.H. Press, B. P. Flannery, * '* S.A. Teukolsky and W.T. Vetterling, Cambridge * '* University Press, 1986" [BIBLI 08]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '******************************************************************** 'PROGRAM TEST_SMOOFT DEFINT I-N NMAX = 2048 F$ = "######.######" DIM ysignal(NMAX) DIM ysave(NMAX) 'open input and output files OPEN "c:\temp\smooth.dat" FOR INPUT AS #1 OPEN "c:\temp\smooth.lst" FOR OUTPUT AS #2 'read number of input ysignal points in input file INPUT #1, ndata 'take nearest power of two IF ndata > 2048 THEN ndata = 2048 IF ndata < 2048 AND ndata > 1023 THEN ndata = 1024 IF ndata < 1024 AND ndata > 511 THEN ndata = 512 IF ndata < 512 AND ndata > 255 THEN ndata = 256 IF ndata < 256 AND ndata > 127 THEN ndata = 128 IF ndata < 128 AND ndata > 63 THEN ndata = 64 IF ndata < 64 AND ndata > 31 THEN ndata = 32 IF ndata < 32 THEN PRINT #2, " Error: number of points too small (<32) !" CLOSE #2 PRINT " Results in file smooth.lst (error)." END END IF 'read ndata couples T(i), Y(i) in input data file FOR i = 1 TO ndata INPUT #1, temp, ysignal(i) IF i = 1 THEN tbegin = temp IF i = ndata THEN tend = temp NEXT i CLOSE #1 FOR i = 1 TO ndata ysave(i) = ysignal(i) NEXT i PTS = 30! GOSUB 1000 'call SMOOFT(ysignal,ndata,30) dt = (tend - tbegin) / (ndata - 1) t = tbegin - dt PRINT #2, " Time Y Smoothed Y " PRINT #2, "----------------------------------------" FOR i = 1 TO ndata t = t + dt PRINT #2, USING F$; t; PRINT #2, USING F$; ysave(i); PRINT #2, USING F$; ysignal(i) NEXT i CLOSE #2 CLS PRINT PRINT " Results in file smooth.lst." PRINT END 'of main program 500 'Function XMAX(a,b) IF a >= b THEN XMAX = a ELSE XMAX = b END IF RETURN 1000 'Subroutine SMOOFT(ysignal,N,PTS) '--------------------------------------------------------------- ' Smooths an array Y of length N, with a window whose full width ' is of order PTS neighboring points, a user supplied value. ' Array Y (here called ysignal) is modified. '--------------------------------------------------------------- m = 2 n = ndata NMIN = n + INT(2 * PTS) 1 IF m < NMIN THEN m = 2 * m GOTO 1 END IF IF m > NMAX THEN PRINT " Sample too big." RETURN END IF CNST = (PTS / m) ^ 2 Y1 = ysignal(1) YN = ysignal(n) RN1 = 1! / (n - 1) FOR j = 1 TO n 'Remove linear trend ysignal(j) = ysignal(j) - RN1 * (Y1 * (n - j) + YN * (j - 1)) NEXT j IF n + 1 <= m THEN FOR j = n + 1 TO m ysignal(j) = 0! NEXT j END IF MO2 = m / 2: isign = 1 GOSUB 2000 'call REALFT(ysignal,MO2,1) Fourier transform ysignal(1) = ysignal(1) / MO2 FAC = 1! FOR j = 1 TO MO2 - 1 K = 2 * j + 1 IF FAC <> 0! THEN 'FAC=XMAX(0.0,(1.0-CNST*J*J)/MO2) a = 0!: b = (1! - CNST * j * j) / MO2: GOSUB 500 ysignal(K) = XMAX * ysignal(K) ysignal(K + 1) = XMAX * ysignal(K + 1) ELSE ysignal(K) = 0! ysignal(K + 1) = 0! END IF NEXT j 'FAC=XMAX(0.0,(1.0-0.25*PTS*PTS)/MO2) Last point a = 0!: b = (1! - .25 * PTS * PTS) / MO2: GOSUB 500 ysignal(2) = XMAX * ysignal(2) isign = -1 GOSUB 2000 'call REALFT(ysignal,MO2,-1) Inverse Fourier transform FOR j = 1 TO n 'Restore linear trend ysignal(j) = RN1 * (Y1 * (n - j) + YN * (j - 1)) + ysignal(j) NEXT j RETURN 1500 'Subroutine four1(ysignal, nn, isign) '-------------------------------------------------------------------------------------------- 'Replaces ysignal(2*nn) by its discrete Fourier transform, if isign is input as 1 or replaces 'ysignal by nn times its inverse discrete Fourier transform, if isign is input as -1. 'ysignal is a complex array of length nn or, equivalently, a real array of length 2*nn. nn 'MUST be an integer power of 2 (this is not checked for here). '-------------------------------------------------------------------------------------------- DIM theta AS DOUBLE 'Double precision for the trigonometric recurrences DIM wi AS DOUBLE DIM wpi AS DOUBLE DIM wpr AS DOUBLE DIM wr AS DOUBLE DIM wtemp AS DOUBLE nn = n / 2 n1 = 2 * nn j = 1 i = 1 WHILE i <= n1 'This is the bit-reversal section of the routine IF j > i THEN tempr = ysignal(j) 'Exchange the two complex numbers tempi = ysignal(j + 1) ysignal(j) = ysignal(i) ysignal(j + 1) = ysignal(i + 1) ysignal(i) = tempr ysignal(i + 1) = tempi END IF m = nn 10 IF (m >= 2 AND j > m) THEN j = j - m m = m / 2 GOTO 10 END IF j = j + m i = i + 2 WEND mmax = 2 'Here begins the Danielson-Lanczos section of the routine 20 IF n1 > mmax THEN 'Outer loop executed log2 nn times istep = 2 * mmax theta = 6.28318530717959# / (isign * mmax) 'Initialize for the trigonometric recurrence wpr = -2# * SIN(.5 * theta) * SIN(.5 * theta) wpi = SIN(theta) wr = 1# wi = 0# m = 1 WHILE m <= mmax 'Here are the two nested inner loops i = m WHILE i <= n1 j = i + mmax 'This is the Danielson-Lanczos formula: tempr = wr * ysignal(j) - wi * ysignal(j + 1) tempi = wr * ysignal(j + 1) + wi * ysignal(j) ysignal(j) = ysignal(i) - tempr ysignal(j + 1) = ysignal(i + 1) - tempi ysignal(i) = ysignal(i) + tempr ysignal(i + 1) = ysignal(i + 1) + tempi i = i + istep WEND wtemp = wr 'Trigonometric recurrence wr = wr * wpr - wi * wpi + wr wi = wi * wpr + wtemp * wpi + wi m = m + 2 WEND mmax = istep GOTO 20 'Not yet done END IF 'All done RETURN 2000 'Subroutine realft(ysignal, n, isign) '-------------------------------------------------------------------------------------------- 'USES four1 'Calculates the Fourier transform of a set of n real-valued data points. Replaces this data '(which is stored in array ysignal(n)) by the positive frequency half of its complex Fourier 'transform. The real-valued first and last components of the complex transform are returned 'as elements ysignal(1) and ysignal(2), respectively. n must be a power of 2. This routine 'also calculates the inverse transform of a complex data array if it is the transform of real 'data (Result in this case must be multiplied by 2/n). '-------------------------------------------------------------------------------------------- DIM PI AS DOUBLE PI = 4# * ATN(1#) n = MO2 theta = PI / (n / 2) 'Initialize the recurrence c1 = .5 IF isign = 1 THEN c2 = -.5 GOSUB 1500 'call four1(ysignal,n/2,1) The forward transform is here ELSE c2 = .5 'Otherwise set up for an inverse transform theta = -theta END IF wpr = -2# * SIN(.5 * theta) * SIN(.5 * theta) wpi = SIN(theta) wr = 1# + wpr wi = wpi n2p3 = n + 3 FOR i = 2 TO n / 4 'Case i=1 done separately below i1 = 2 * i - 1 i2 = i1 + 1 i3 = n2p3 - i2 i4 = i3 + 1 wrs = wr wis = wi h1r = c1 * (ysignal(i1) + ysignal(i3))'The two separate transforms are separated out of data h1i = c1 * (ysignal(i2) - ysignal(i4)) h2r = -c2 * (ysignal(i2) + ysignal(i4)) h2i = c2 * (ysignal(i1) - ysignal(i3)) ysignal(i1) = h1r + wrs * h2r - wis * h2i'Here they are recombined to form the true transform 'of the original real data. ysignal(i2) = h1i + wrs * h2i + wis * h2r ysignal(i3) = h1r - wrs * h2r + wis * h2i ysignal(i4) = -h1i + wrs * h2i + wis * h2r wtemp = wr 'The recurrence wr = wr * wpr - wi * wpi + wr wi = wi * wpr + wtemp * wpi + wi NEXT i IF isign = 1 THEN h1r = ysignal(1) ysignal(1) = h1r + ysignal(2) ysignal(2) = h1r - ysignal(2) 'Squeeze the first and last data together to get ELSE 'them all within the original array h1r = ysignal(1) ysignal(1) = c1 * (h1r + ysignal(2)) ysignal(2) = c1 * (h1r - ysignal(2)) GOSUB 1500 'call four1(ysignal,n/2,-1) This is the inverse transform for the case isign=-1 END IF RETURN 'end of file smooth.bas