10   '
20   '            FOURIER SYNTHESIS OF WAVEFORMS
30   '
40   '--------------------- 1/5/88 -------------------------
50   '
60   '  This program does a Fourier synthesis of four
     '  functions.  It adds some number of sine waves, each
     '  multiplied by the appropriate coefficient, and makes
     '  a plot of the resultant wave form.  The four
     '  functions programmed are square wave, triangle,
     '  spike, and face profile.  It also calculates and
     '  plots the average deviation for 1, 2, 3, ... , N
     '  waves if not zero and if the wave is not the
     '  profile.
130  '
140  '           --> LIST OF PROCEDURES <--
150  '
160        GOSUB 340   'INITIALIZATION 1
170        GOSUB 540   'READ IN VALUES FOR THIS RUN
180        GOSUB 760   'INITIALIZATION 2
190        GOSUB 880   'CALC. SERIES APPROX. FOR EACH X
200        GOSUB 1380  'PRINT THE PARAMETERS FOR THIS RUN
210        GOSUB 1510  'PLOT APPROX. FUNC. VS. X
220  '
230  PRINT
240  INPUT "PLOT AV. DEV. vs NO. TERMS (Y/N)? ", ANS$
250  IF ANS$ = "y" THEN ANS$ = "Y"
260  IF ANS$ = "Y" THEN GOSUB 2070  'PLOT AV. DEV. vs NO._
     TERMS
270  '
280  PRINT
290  INPUT "WANT TO RUN ANOTHER (Y/N)? ", AGAIN$
300  IF AGAIN$ = "Y" OR AGAIN$ = "y" THEN CLOSE #1: GOTO 170
310  '
320  END
330  '------------------------------------------------------
340  'INITIALIZATION 1
350  '
360  DIM Z$(101), DS(100), A(48), B(48), XS(101), FUNCS(101)
370  PI = 3.14159265#
380  READ BLANK$, VLN$, ASTER$, DOT$
390  DATA  " ","|","*","."
400  '
410  GOSUB 3020  'STORE THE COEFFICIENTS FOR THE PROFILE
420  '
430  RM = 63        ' "RM" is the (odd) graph right margin;_
     63 or 101
440  LM = 17        ' "LM" is the graph left margin
450  ZM = INT((RM - 1) / 2) + 1  
460  '"ZM" is the M value of the vert. center line
470  ' -- Initialize Z$ to blanks:
480  FOR I = 1 TO RM
490    Z$(I) = BLANK$
500    NEXT I
510  '
520  RETURN
530  '------------------------------------------------------
540  'READ IN VALUES FOR THIS RUN
550  '
560  CLS
570  PRINT "INPUT VALUES FOR THE FOLLOWING:"
580  PRINT "  NUMBER OF WAVES (N)"
590  PRINT "  WAVE TYPE (FLAG)"
600  ' -- FLAG: 0 FOR SQUARE WAVE
610  '          1 FOR TRIANGLE
620  '          2 FOR SPIKE
630  '          3 FOR FACE PROFILE
640  '  
650  PRINT "  NO. OF X POINTS (ODD)"
660  PRINT
670  PRINT "SEPARATE THE FOUR VALUES WITH COMMAS,"
680  PRINT "THEN PRESS RETURN"
690  PRINT
700  PRINT
710  PRINT
720  INPUT N, FLAG, NL
730  '
740  RETURN
750  '------------------------------------------------------
760  'INITIALIZATION 2
770  '
780  N = INT(N)
790  NL = INT(NL)           
          ' "NL" is the (odd) # of vert. graph lines
800  VZ = INT((NL - 1) / 2) + 1         
          ' "VZ" is the line # where X=0
810  SCF = 1.5           
          '"SCF" is scale factor for plotting FUNC vs X
820  IF FLAG = 2 THEN SCF = 1.25
830  IF FLAG = 3 THEN SCF = 1
840  '
850  DX = 2 * PI / (NL - 1)
860  RETURN
870  '------------------------------------------------------
880  'CALC. APPROX. FUNC AT NL X-POINTS
890  '
900  X = -PI
910  FOR J = 1 TO NL
920    '
930    ' -- FUNC (the sum of the N waves) is initialized to_
       '    zero except for the spike wave.
950    '
960    FUNC = 0
970    IF FLAG = 2 THEN FUNC = 1
980    '
990    ' -- Begin loop over # of waves (terms) to be added
       '    in the sum.
1000   '
1010   FOR K = 1 TO N
1020     '
1030     XK = K
1040     '
1050     ' -- DS(K) accumulates average deviation of the
         '    approximate value (for K waves) from exact
         '    value, summed over all X.
1070     '
1080     IF J = 1 THEN DS(K) = 0
1090     '
1100     ' --   AK and BK are the Fourier coefficients
1110     '
1120     AK = 0
1130     BK = 0
1140     IF FLAG = 0 THEN GOSUB 2590  'SQUARE WAVE COEFF'S
1150     IF FLAG = 1 THEN GOSUB 2740  'TRIANGLE WAVE COEFF'S
1160     IF FLAG = 2 THEN GOSUB 2830  'SPIKE WAVE COEFF'S
1170     IF FLAG = 3 THEN GOSUB 2900  'FACE PROFILE COEFF'S
1180     '
1190     ' -- Add this wave to the sum.
1200     '
1210     FUNC = FUNC + AK * COS(XK * X)
1220     FUNC = FUNC + BK * SIN(XK * X)
1230     '
1240     IF FLAG <> 2 THEN DS(K) = DS(K) + ABS(F - FUNC)/NL
1250     IF FLAG = 2 AND J <> VZ THEN DS(K) = DS(K) + 
1260               ABS(F - FUNC) / NL
1270     NEXT K
1280   '
1290   XS(J) = X
1300   FUNCS(J) = FUNC
1310   EN = N
1320   IF FLAG = 2 THEN FUNCS(J) = FUNC / (EN + 1)
1330   X = X + DX
1340   NEXT J
1350 '
1360 RETURN
1370 '-------------------------------------------------------
1380 'PRINT OUT VALUES FOR THIS RUN
1390 '
1400 INPUT "PRINTER (Y/N)? ", ANS$
1410 IF ANS$ = "y" OR ANS$ = "Y" THEN
         OPEN "LPT1:" FOR OUTPUT AS #1
1420 IF ANS$ <> "y" AND ANS$ <> "Y" THEN
         OPEN "SCRN:" FOR OUTPUT AS #1
1430 '
1440 PRINT #1, "NUMBER OF WAVES: "; N
1450 PRINT #1, "WAVE TYPE:  "; FLAG
1460 PRINT #1, "NO. OF X VALUES:"; NL
1470 PRINT #1,
1480 '
1490 RETURN
1500 '------------------------------------------------------
1510 'PLOT APPROX. FUNC. VS. X
1520 '
1530 PRINT #1, "   X      F(X)  ";
1540 FOR I = 1 TO RM
1550   PRINT #1, DOT$;
1560   NEXT I
1570 PRINT #1,
1580 '
1590 FOR J = 1 TO NL
1600   X = XS(J)
1610   FUNC = FUNCS(J)
1620   ' -- Compute an integer in the range M=1,...,RM to
       '    correspond to the value of the resultant wave
       '    amplitude stored in FUNC. Note that when FUNC=0
       '    we have M in the center of the plot.
1660   ADD = .00001
1670   IF FUNC < 0 THEN ADD = .49999
1680   IF FUNC > 0 THEN ADD = .5
1690   M = 1 + INT((ZM - 1) * (FUNC / SCF + 1) + ADD)
1700   IF M < 1 OR M > RM THEN M = 0
1710   '
1720   ' -- If not profile, mark the X-axis
1730   '
1740   IF FLAG <> 3 THEN LET Z$(ZM) = VLN$
1750   Z$(1) = DOT$
1760   Z$(RM) = DOT$
1770   '
1780   ' -- Store an asterisk in Z$(M)
1790   '
1800   Z$(M) = ASTER$
1810   PRINT #1, USING "##.####";_
           INT((X + .00005) *10^4)/10^4;
1820   PRINT #1, " ";
1830   PRINT #1, USING "##.####";_
           INT((FUNC + .00005) *10^4)/10^4;
1840   PRINT #1, " ";
1850   FOR I = 1 TO RM
1860     PRINT #1, Z$(I);
1870     NEXT I
1880   PRINT #1,
1890   '
1900   ' -- Erase the asterisk before moving to the next
       '    point.
1910   '
1920   Z$(M) = BLANK$
1930   '
1940   NEXT J
1950 '
1960 ' -- Put a line of dots across the bottom of the plot,
1970 '    then a blank line.
1980 '
1990 PRINT #1, TAB(LM);
2000 FOR I = 1 TO RM
2010   PRINT #1, DOT$;
2020   NEXT I
2030 PRINT #1,
2040 '
2050 RETURN
2060 '------------------------------------------------------
2070 'PLOT AV. DEV. vs. NO. TERMS
2080 '
2090 FOR NP = 1 TO 10
2100   PRINT #1,
2110   NEXT NP
2120 Z$(ZM) = BLANK$
2130 PRINT #1, "AVERAGE DEVIATION VS. NUMBER OF TERMS"
2140 PRINT #1,
2150 PRINT #1, "      N   DEV   ";
2160 FOR I = 1 TO RM
2170   PRINT #1, DOT$;
2180   NEXT I
2190 PRINT #1,
2200 '
2210 ' FIND AV. DS FOR SCALE FACTOR
2220 SUM = 0
2230 IF FLAG <> 2 THEN 2280
2240 FOR K = 1 TO N
2250   XK = K
2260   DS(K) = DS(K) / (XK + 1)
2270   NEXT K
2280 FOR K = 2 TO N
2290   SUM = SUM + DS(K)
2300   NEXT K
2310 SCF = (ZM - 1) * (N - 1) / SUM
2320 '
2330 ' BEGIN PLOT LOOP
2340 FOR K = 2 TO N
2350   M = INT(SCF * DS(K) + .49999) + 1
2360   IF M < 1 OR M > RM THEN M = 0
2370   Z$(1) = DOT$
2380   Z$(RM) = DOT$
2390   Z$(M) = ASTER$
2400   PRINT #1, USING "#######"; K;
2410   PRINT #1, " ";
2420   PRINT #1, USING "##.####";_
           INT((DS(K) + .00005) *10^4)/10^4;
2430   PRINT #1, " ";
2440   FOR I = 1 TO RM
2450     PRINT #1, Z$(I);
2460     NEXT I
2470   PRINT #1,
2480   Z$(M) = BLANK$
2490   NEXT K
2500 '
2510 PRINT #1, "                ";
2520 FOR I = 1 TO RM
2530   PRINT #1, DOT$;
2540   NEXT I
2550 PRINT #1,
2560 '
2570 RETURN
2580 '------------------------------------------------------
2590 'SQUARE WAVE
2600 '
2610 ' -- F is the theoretical function (in this case square
     '    wave).
2620 '
2630 F = 0
2640 IF X <> 0 THEN LET F = ABS(X) / X
2650 '
2660 ' -- The factor C makes AK zero for even values of K.
2670 '
2680 C = K - 2 * INT(K / 2)
2690 BK = 0
2700 IF XK <> 0 THEN BK = 4 * C / (PI * XK)
2710 '
2720 RETURN
2730 '------------------------------------------------------
2740 'TRIANGULAR WAVE
2750 '
2760 F = 1 - 2 * ABS(X) / PI
2770 C = K - 2 * INT(K / 2)
2780 AK = 0
2790 IF XK <> 0 THEN AK = 8 * C / (PI * XK) ^ 2
2800 '
2810 RETURN
2820 '------------------------------------------------------
2830 'SPIKE FUNCTION
2840 '
2850 F = 0
2860 AK = 1
2870 '
2880 RETURN
2890 '------------------------------------------------------
2900 'FACE PROFILE
2910 '
2920 ' -- Only the first 48 coefficients are available for
     '    the profile.
2940 AK = 0
2950 BK = 0
2960 IF K > 48 THEN RETURN
2970 AK = A(K)
2980 BK = B(K)
2990 '
3000 RETURN
3010 '------------------------------------------------------
3020 'STORE THE COEFFICIENTS FOR THE FACE PROFILE
3030 '
3040 FOR I = 1 TO 48
3050   READ A(I)
3060   NEXT I
3070 '
3080 DATA  .1322, .0138, .0816,-.0360, .0122,-.0379
3090 DATA -.0046,-.0147, .0089, .0036, .0084, .0050
3100 DATA -.0027,-.0015,-.0055, .0013,-.0042, .0050
3110 DATA -.0017, .0030,-.0041, .0010,-.0047,-.0002
3120 DATA -.0019, .0020, .0001, .0017, .0002,-.0011
3130 DATA -.0016,-.0021,-.0001,-.0008, .0023, .0002
3140 DATA  .0015,-.0016, .0003,-.0016, .0001,-.0002
3150 DATA  .0013, .0004, .0007,-.0003,-.0006,-.0009
3160 FOR I = 1 TO 48
3170   READ B(I)
3180   NEXT I
3190 '
3200 DATA -.1534,-.0468, .0777,-.0440, .0593,-.0544
3210 DATA  .0387,-.0457, .0348,-.0322, .0381,-.0229
3220 DATA  .0306,-.0210, .0205,-.0246, .0133,-.0202
3230 DATA  .0127,-.0130, .0142,-.0103, .0099,-.0121
3240 DATA  .0066,-.0117, .0075,-.0081, .0088,-.0081
3250 DATA  .0079,-.0090, .0065,-.0093, .0081,-.0081
3260 DATA  .0093,-.0074, .0089,-.0092, .0076,-.0098
3270 DATA  .0074,-.0088, .0091,-.0077, .0084,-.0088
3280 '
3290 RETURN
3300 '---------------- END OF PROGRAM LINES ----------------