'===========================================================================
' Subject: CALCULATE FOURIER COEFFICIENTS    Date: 10-04-99 (02:50)       
' Author:  David Karlov                      Code: QB, QBasic, PDS        
' Origin:  dkarlov@acs.itd.uts.edu.au      Packet: ALGOR.ABC
'===========================================================================
' This program can calculate the Fourier coefficients of a periodic function.
' If you don't know what that means, a maths textbook can help..
' But it isn't necessary to understand it to use this program.

' Here is a very brief description of what all this Fourier stuff is:
' Every function can be thought of as a sum of an infinite number of
' sinusoidal functions.  Each sinusoid has a frequency which is an integer
' multiple of the "fundamental frequency", or the frequency of the original
' function (1 over the period).  These sinusoids are added together, and by
' changing the amplitude (strength) of each sinusoid, any periodic function
' can be made.
' This program finds what amplitudes are necessary to make a few interesting
' functions.  You can look at a table or graph of those values, or watch the
' function being built step-by-step out of sinusoids, or view an animated
' screen where each sinusoid is a rotating line attached on to the end of
' the previous rotating line, each rotating at different speeds (the speed
' represents the frequency of the sinusoid).  All of these line segments
' rotating together make up the original function ..approximately ;)
'
' NOTE -- if you want to add your own functions to do stuff with, look in
' the QB function F1(), add the code in there, and alter the constant
' NumFuncs.  You have to specify the function mathematically so QB can
' understand it, and also state the function's period.
'
' ..Written by David Karlov
' Feel free to contact me:  dkarlov@rocketmail.com


TYPE GraphData
   MinX AS DOUBLE
   MaxX AS DOUBLE
   MinY AS DOUBLE
   MaxY AS DOUBLE
   XRange AS DOUBLE
   YRange AS DOUBLE
END TYPE

DECLARE FUNCTION Mod2# (x AS DOUBLE, t AS DOUBLE)
DECLARE FUNCTION ArcTan# (Numer AS DOUBLE, Denom AS DOUBLE)
DECLARE FUNCTION Rect# (x AS DOUBLE)
DECLARE FUNCTION Sinc# (x AS DOUBLE)
DECLARE FUNCTION SquareWave# (x AS DOUBLE, DC AS DOUBLE)

DECLARE FUNCTION IntegrateF1# (x0 AS DOUBLE, x1 AS DOUBLE, Condition AS INTEGER, n AS INTEGER)

DECLARE SUB ShowPage (PageNum AS INTEGER)
DECLARE SUB GraphF1 (Params AS GraphData)
DECLARE SUB GraphAxes (Params AS GraphData)
DECLARE SUB GraphFourier (Params AS GraphData, Num AS INTEGER)
DECLARE SUB ShowFourierCoeffs ()
DECLARE SUB GraphFourierCoeffs ()
DECLARE SUB SubGraph (d() AS DOUBLE, x0 AS INTEGER, y0 AS INTEGER, Label AS STRING)
DECLARE SUB DrawWhirly ()
DECLARE SUB CalculateCoeffs ()
DECLARE SUB Options ()
DECLARE FUNCTION Trim$ (a AS DOUBLE)

DECLARE FUNCTION F1# (x AS DOUBLE, Condition AS INTEGER, n AS INTEGER)

CONST IntegrateResolution = 500  ' The higher this is, the more accurate
                                 ' (but slower also).
CONST ScreenMode = 9
CONST XSize = 640
CONST YSize = 350
CONST ExtraYSpace = .2
CONST Pi = 3.141592653589793#
CONST Tiny = .00000000000001#
CONST NUMCOEFFS = 70          ' Number of Fourier coefficients to calculate.

CONST NumFuncs = 12
DIM SHARED FF AS INTEGER ' the current function number. (1..NumFuncts)

DIM SHARED T0 AS DOUBLE  ' period of function


DIM SHARED FourierCoeffs(1 TO 2, 0 TO NUMCOEFFS) AS DOUBLE
DIM SHARED Func$
DIM SHARED Params AS GraphData

DIM SHARED Page AS INTEGER
Page = 0: ShowPage Page

FF = 1   ' Start on function 1

CalculateCoeffs

Options

' Inverse TAN function. (returns result in radians)
FUNCTION ArcTan# (Numer AS DOUBLE, Denom AS DOUBLE)

   DIM i AS DOUBLE

   IF Numer = 0 THEN
      IF Denom >= 0 THEN ArcTan = 0 ELSE ArcTan = Pi
   ELSEIF Denom = 0 THEN
      IF Numer > 0 THEN ArcTan = Pi / 2 ELSE ArcTan = -Pi / 2
   ELSE
      i = ATN(Numer / Denom)
      IF SGN(Numer) * SGN(Denom) = -1 THEN i = i + Pi
      ArcTan = i
   END IF

END FUNCTION

' Calculates the Fourier coefficients for the current function.
SUB CalculateCoeffs

   DIM n AS INTEGER
   n = F1(0, 0, 0) ' dummy call to initialize T0.

   FOR n = 0 TO NUMCOEFFS
      IF n = 0 THEN
         FourierCoeffs(1, n) = IntegrateF1(-T0 / 2, T0 / 2, 0, 0) / T0
      ELSE
         FourierCoeffs(1, n) = IntegrateF1(-T0 / 2, T0 / 2, 1, n) * 2 / T0
      END IF
      FourierCoeffs(2, n) = IntegrateF1(-T0 / 2, T0 / 2, 2, n) * 2 / T0
   NEXT n

END SUB

' Draws the rotating phasors on the screen.
' The image is drawn as fast as your computer can do it, with no delay.
SUB DrawWhirly
   DIM Im1(30000) AS INTEGER
 
   SCREEN 9, 0, 0, 0: CLS
   SCREEN 9, 0, 1, 1: CLS

   AngularRes = 100
   DetailRes = 70
   Scale1 = 175

   EX1 = 0: Beg = 0: Oldy = -1
   GET (1, 0)-(290, 349), Im1(0)

   DO
      FOR Ang = 1 TO AngularRes
         Page = 1 - Page: 'ShowPage Page
         SCREEN ScreenMode, , Page, Page
         PCOPY Page, 1 - Page

         IF Beg = 1 THEN DRAW "C1A0L640"
        
         FOR y = 0 TO 349
            IF POINT(0, y) = 1 THEN EXIT FOR
         NEXT y

         Page = 1 - Page: ShowPage Page
         LINE (290, 0)-(290, 349), 0
         PUT (0, 0), Im1(0), PSET
         IF y < 350 THEN
            IF Oldy > -1 THEN
               LINE (289, Oldy)-(290, y), 15
            ELSE
               PSET (290, y), 15
            END IF
            Oldy = y
         END IF
         Beg = 1

         GET (1, 0)-(290, 349), Im1(0)

         PRESET (320, 175)
         DRAW "C14"
         FOR t = 1 TO DetailRes
            a = FourierCoeffs(1, t)
            b = FourierCoeffs(2, t)
            c = (-ArcTan((b), (a)) * 180 / Pi + t * Ang / AngularRes * 360) MOD 360
            DRAW "TA=" + VARPTR$(c)
            c = INT(Scale1 * SQR(a * a + b * b))
            DRAW "U=" + VARPTR$(c)
         NEXT t

         DRAW "C15U1L1D2R2U2L1"
         IF INKEY$ <> "" THEN EX1 = 1: EXIT FOR
      NEXT Ang
   LOOP UNTIL EX1 = 1

END SUB

' This function is the actual function we are finding the Fourier
' coefficients for.
' Depending on the value of FF (a shared variable), one of a number of
' different functions are selected, each of which is a function of x.
'
' The input Condition is used to help find the Fourier coefficients..
' 0 means the actual function value is returned
' 1 means the value is multiplied by COS(2 * Pi * n / T0 * x)
' 2 means the value is multiplied by SIN(2 * Pi * n / T0 * x)
'
' (T0 is the period of the function, it is a shared variable).
'
FUNCTION F1# (x AS DOUBLE, Condition AS INTEGER, n AS INTEGER)

   SELECT CASE FF
   CASE 1
      T0 = 4   ' Period
      i# = SquareWave((x - 1) / 4, 2 / 4) - .5  ' Actual function
      Func$ = "g(t) =  Sigma(n = -oo -> n = oo) Rect((t - 1 - 4n) / 2)- 1/2"
      ' Func$ can be used to hold a description of the function.
   CASE 2
      T0 = 4
      IF x = 0 THEN
         i# = .5
      ELSE
         i# = SquareWave((x - 3) / 4, .5) - .5
      i# = i# + .5 * SquareWave((x - 2.5) / 4, .25) * (Mod2(x, T0) - 3)
         i# = i# + .5 * SquareWave((x - .5) / 4, .25) * (1 - Mod2(x, T0))
      END IF
      Func$ = ""
   CASE 3
      T0 = 4
      i# = SquareWave(x / 4, 2 / 4) * COS(4 * Pi * x)
      Func$ = "g(t) = Sigma(n = -oo -> n = oo) Rect(t - 4n) / 2) * Cos(4 Pi t)"
   CASE 4
      T0 = 2 * Pi
      i# = SquareWave((x - Pi / 2) / (2 * Pi), .5) * EXP(-Mod2(x, T0)) / 2
      i# = i# + (SquareWave((x - Pi / 2) / (2 * Pi), .5) - 1) * EXP(-Mod2(x, T0 / 2)) / 2
      Func$ = ""
   CASE 5
      T0 = 2
      IF x = 0 THEN
         i# = 0
      ELSE
         i# = SquareWave((x - T0 / 4) / T0, .5) * Mod2(x, T0)
         i# = i# + SquareWave((x - 3 * T0 / 4) / T0, .5) * (2 - Mod2(x, T0))
      END IF
      Func$ = ""
   CASE 6
      T0 = 2
      IF x = 0 THEN
         i# = 0
      ELSE
         i# = SquareWave((x - T0 / 4) / T0, .5) * Mod2(x, T0)
         i# = i# + SquareWave((x - 3 * T0 / 4) / T0, .5) * (2 - Mod2(x, T0))
      END IF
      i# = i# * (SIN(2 * Pi * x))
      Func$ = ""
   CASE 7
      T0 = 4
      i# = .5 * SIN(EXP(.8 * Mod2(x, T0)))
      Func$ = "w(t) = Sigma (n = -oo -> n = oo) Rect((t - 4n) / 4) * Sin(.8 * e ^ (t - 4n))"
   CASE 8
      T0 = 1 * 20
      i# = COS(2 * Pi * x / 20) * SquareWave(x * 3 * 2 * Pi / 20, .1)
      Func$ = ""
   CASE 9
      T0 = 20
      i# = COS(2 * Pi * x / 20) * (COS(2 * Pi * 33 * x / 20)) ^ 2
      Func$ = ""
   CASE 10
      T0 = 2
      i# = SQR(1 - (Mod2(x, T0) - 1) ^ 2)
      Func$ = ""
   CASE 11
      T0 = 2
      i# = .2 * (SIN(Pi * x) + SIN(3 * Pi * x) + SIN(5 * Pi * x) + SIN(7 * Pi * x))
      Func$ = "w(t) = (Sin(Pi t) + Sin(3 Pi t) + Sin(5 Pi t) + Sin(7 Pi t)) / 5"
   CASE 12
      T0 = 2
      i# = (SquareWave((x - 1) / 2, .5) * COS(Pi * x) + ABS(SquareWave(x / 2, .5) * SIN(2 * Pi * x)) / 2)
      Func$ = ""
   END SELECT

   SELECT CASE Condition
      CASE 0   ' normal
      CASE 1   ' projected onto COS
         i# = i# * COS(2 * Pi * n / T0 * x)
      CASE 2   ' projected onto SIN
         i# = i# * SIN(2 * Pi * n / T0 * x)
   END SELECT
  
   F1 = i#

END FUNCTION

' Draw the axes for the original function (in options menu).
SUB GraphAxes (Params AS GraphData)

   ' uses F1 as a guide for the size of the axes.

   DIM Y1(0 TO XSize - 1) AS DOUBLE

   IF Params.XRange = 0 THEN
      Params.MinX = -10.5: Params.MaxX = 10.5
      Params.XRange = Params.MaxX - Params.MinX
      Params.MinY = 30000: Params.MaxY = -30000

      FOR i = 0 TO XSize - 1
         Y1(i) = F1((i / XSize) * Params.XRange + Params.MinX, 0, 0)
         IF Y1(i) > Params.MaxY THEN Params.MaxY = Y1(i)
         IF Y1(i) < Params.MinY THEN Params.MinY = Y1(i)
      NEXT i
 
      Params.YRange = Params.MaxY - Params.MinY
      Params.MaxY = Params.MaxY + Params.YRange * ExtraYSpace
      Params.MinY = Params.MinY - Params.YRange * ExtraYSpace

      Params.YRange = Params.MaxY - Params.MinY
   END IF

   i = -XSize * Params.MinX / Params.XRange
   LINE (i, 0)-(i, YSize - 1), 1

   i = YSize * (1 + Params.MinY / Params.YRange)
   LINE (0, i)-(XSize - 1, i), 1

   XStep = 10 ^ INT(LOG(Params.XRange) / 3.5)
   FOR i = 0 TO Params.MinX STEP -XStep
      LINE ((i - Params.MinX) / Params.XRange * XSize, YSize * (1 + Params.MinY / Params.YRange) - 5)-((i - Params.MinX) / Params.XRange * XSize, YSize * (1 + Params.MinY / Params.YRange) + 5), 1
   NEXT i
   FOR i = 0 TO Params.MaxX STEP XStep
      LINE ((i - Params.MinX) / Params.XRange * XSize, YSize * (1 + Params.MinY / Params.YRange) - 5)-((i - Params.MinX) / Params.XRange * XSize, YSize * (1 + Params.MinY / Params.YRange) + 5), 1
   NEXT i

   YStep = 10 ^ INT(LOG(Params.YRange) / 3.5)
   FOR i = 0 TO Params.MinY STEP -YStep
      LINE (-XSize * Params.MinX / Params.XRange - 5, YSize - 1 - (i - Params.MinY) / Params.YRange * YSize)-(-XSize * Params.MinX / Params.XRange + 5, YSize - 1 - (i - Params.MinY) / Params.YRange * YSize), 1
   NEXT i
   FOR i = 0 TO Params.MaxY STEP YStep
      LINE (-XSize * Params.MinX / Params.XRange - 5, YSize - 1 - (i - Params.MinY) / Params.YRange * YSize)-(-XSize * Params.MinX / Params.XRange + 5, YSize - 1 - (i - Params.MinY) / Params.YRange * YSize), 1
   NEXT i


END SUB

' Draw the actual function (in options menu).
SUB GraphF1 (Params AS GraphData)

   DIM Y1(0 TO XSize - 1) AS DOUBLE

   FOR i = 0 TO XSize - 1
      Y1(i) = F1((i / XSize) * Params.XRange + Params.MinX, 0, 0)
   NEXT i

   PSET (0, (Params.MaxY - Y1(0)) / Params.YRange * YSize), 12
   FOR i = 1 TO XSize - 1
      LINE -(i, (Params.MaxY - Y1(i)) / Params.YRange * YSize), 12
   NEXT i

END SUB

' Draw the Fourier transform of the function, using Num sinsoids.
SUB GraphFourier (Params AS GraphData, Num AS INTEGER)

   DIM Y1(0 TO XSize - 1) AS DOUBLE

   FOR i = 0 TO XSize - 1
      Y1(i) = 0
      FOR j = 0 TO Num
         Y1(i) = Y1(i) + FourierCoeffs(1, j) * COS(2 * Pi * j / T0 * ((i / XSize) * Params.XRange + Params.MinX))
         Y1(i) = Y1(i) + FourierCoeffs(2, j) * SIN(2 * Pi * j / T0 * ((i / XSize) * Params.XRange + Params.MinX))
      NEXT j
   NEXT i

   PSET (0, (Params.MaxY - Y1(0)) / Params.YRange * YSize), 14
   FOR i = 1 TO XSize - 1
      LINE -(i, (Params.MaxY - Y1(i)) / Params.YRange * YSize), 14
   NEXT i


END SUB

' Draw 4 graphs, showing the spectra of the function.
' Two of them are a(n) and b(n), these are the coefficients of the COS and
' SIN terms of the Fourier series.  A(n) and Phi(n) form the compact Fourier
' series, A(n) being the magnitude and Phi(n) the phase for each n.
'
' These graphs are another way to show the tables in ShowFourierCoeffs().
SUB GraphFourierCoeffs

   SCREEN 12: CLS

   DIM d(0 TO 50) AS DOUBLE   ' data array

   ' do 4 graphs:    a(n) vs n
   '                 b(n) vs n
   '                 A(n) vs n
   '                 Phi(n) vs n

   ' GRAPH: a(n) vs n in top left corner
   FOR i = 0 TO 50: d(i) = FourierCoeffs(1, i): NEXT i
   SubGraph d(), 60, 20, "a(n)"

   ' GRAPH: b(n) vs n in bottom left corner
   FOR i = 0 TO 50: d(i) = FourierCoeffs(2, i): NEXT i
   SubGraph d(), 60, 260, "b(n)"

   ' GRAPH: A(n) vs n in top right corner
   FOR i = 0 TO 50
      d(i) = SQR(FourierCoeffs(1, i) ^ 2 + FourierCoeffs(2, i) ^ 2)
   NEXT i
   SubGraph d(), 380, 20, "A(n)"

   ' GRAPH: Phi(n) vs n in bottom right corner
   FOR i = 0 TO 50
      j1 = FourierCoeffs(2, i): IF ABS(j1) < .000001 THEN j1 = 0
      j2 = FourierCoeffs(1, i): IF ABS(j2) < .000001 THEN j2 = 0
      d(i) = -ArcTan((j1), (j2))
      IF d(i) < -Pi THEN d(i) = d(i) + 2 * Pi
      IF d(i) > Pi THEN d(i) = d(i) - 2 * Pi
   NEXT i
   SubGraph d(), 380, 260, "Phi(n)"

   LINE (319, 0)-(319, 463), 8
   LINE (0, 239)-(639, 239), 8
   LINE (0, 0)-(639, 463), 8, B
   LOCATE 30, 2: COLOR 7: PRINT Func$;
  
   WHILE INKEY$ = "": WEND
  
END SUB

' Performs a crude numerical integration of the function specified by FF,
' between x0 and x1.
' The values Condition and n are used by F1() to simplify finding Fourier
' coefficients.
FUNCTION IntegrateF1# (x0 AS DOUBLE, x1 AS DOUBLE, Condition AS INTEGER, n AS INTEGER)

   DIM XRange AS DOUBLE, Total AS DOUBLE
   DIM i AS INTEGER

   XRange = x1 - x0
   FOR i = 0 TO IntegrateResolution - 1
      Total = Total + F1((i / IntegrateResolution) * XRange + x0, Condition, n)
   NEXT i

   IntegrateF1 = Total * XRange / IntegrateResolution

END FUNCTION

' Mod2 is a lot like the built-in MOD function, except that with Mod2,
' the values do not have to be integers.
' e.g.  Mod2(4.6, 1.1) will return 0.2
FUNCTION Mod2# (x AS DOUBLE, t AS DOUBLE)
  
   DIM NewX AS DOUBLE

   NewX = x - t * INT(x / t)
   DO UNTIL NewX < 0: NewX = NewX - t: LOOP
   DO UNTIL NewX >= 0: NewX = NewX + t: LOOP
   Mod2 = NewX

END FUNCTION

' Draw a graph of the function and ask the user what to do next.
SUB Options

   DIM n AS INTEGER

   DO
      GraphAxes Params
      GraphF1 Params
      COLOR 15
      LOCATE 1, 40 - LEN(Func$) / 2: PRINT Func$;
      COLOR 7
      LOCATE 20, 1: PRINT " 1. Select a different function";
      LOCATE 21, 1: PRINT " 2. Show table of Fourier coefficients";
      LOCATE 22, 1: PRINT " 3. Show spectra of Fourier coefficients";
      LOCATE 23, 1: PRINT " 4. Build function from sinusoids";
      LOCATE 24, 1: PRINT " 5. Watch rotating phasors";
      LOCATE 25, 1: PRINT " 6. Quit";
      Page = 1 - Page: ShowPage Page
      DO: i$ = INKEY$: LOOP WHILE i$ = ""
      SELECT CASE i$
         CASE "1"
            SCREEN 9, , Page, Page
            DO
               FOR y = 20 TO 25
                  LOCATE y, 1: PRINT SPACE$(80);
               NEXT y
               LOCATE 22, 1: PRINT " There are"; NumFuncs; "defined functions."
               PRINT " Enter a number from 1 to"; NumFuncs; "-> ";
               INPUT F$
               v = VAL(F$)
            LOOP WHILE v < 1 OR v > NumFuncs OR v <> INT(v)
            FF = v
            Params.XRange = 0    ' Force recalculation of graph parameters
            PRINT : PRINT " Please wait, calculating Fourier coefficients....";
            CalculateCoeffs
            Page = 1 - Page: ShowPage Page
         CASE "2"
            ShowFourierCoeffs
            Page = 1 - Page: ShowPage Page
         CASE "3"
            GraphFourierCoeffs
            Page = 1 - Page: ShowPage Page
         CASE "4"
            CLS
            Page = 1 - Page: ShowPage Page
            CLS
            FOR n = 0 TO NUMCOEFFS
               GraphAxes Params
               GraphF1 Params
               GraphFourier Params, n
               LOCATE 1, 1: PRINT "n ="; n
               IF n > 0 THEN
                  DO: i$ = INKEY$: LOOP UNTIL i$ <> ""
                  IF i$ = CHR$(27) THEN Page = 1 - Page: ShowPage Page: EXIT FOR
               END IF
               Page = 1 - Page: ShowPage Page
            NEXT n
         CASE "5"
            DrawWhirly
            Page = 1 - Page: ShowPage Page
         CASE "6"
            SCREEN 0, 0, 0, 0: CLS : END
      END SELECT
   LOOP


END SUB

' A basic building block to make more complicated functions out of.
' Rect(x) = 1 when x is between -0.5 and 0.5, and 0 elsewhere.
'
' Using the form Rect((x - x0) / A), this returns 1 when the difference
' between x and x0 is less than A/2, and 0 otherwise.
'
FUNCTION Rect# (x AS DOUBLE)

   IF ABS(x) > .5 THEN
      Rect = 0
   ELSE
      Rect = 1
   END IF

END FUNCTION

' Draw 4 tables, showing the spectra of the function.
' Two of them are a(n) and b(n), these are the coefficients of the COS and
' SIN terms of the Fourier series.  A(n) and Phi(n) form the compact Fourier
' series, A(n) being the magnitude and Phi(n) the phase for each n.
'
' These tables show the values used for the graphs in GraphFourierCoeffs().
SUB ShowFourierCoeffs

   SCREEN 0, 0, 0, 0: WIDTH 80, 50: CLS

   DIM a AS SINGLE, b AS SINGLE
  
   LOCATE 1, 1
   FOR n = 0 TO 49
      LOCATE n + 1, 1
      a = FourierCoeffs(1, n)
      IF ABS(a) < Tiny THEN a = 0
      b = FourierCoeffs(2, n)
      IF ABS(b) < Tiny THEN b = 0
      Num1$ = Trim$((n))
      COLOR 15: LOCATE , 1: PRINT "a";
      COLOR 7: PRINT Num1$; " = ";
      COLOR 1: PRINT USING "\    \"; Trim$((a));
      COLOR 15: LOCATE , 18: PRINT "b";
      COLOR 7: PRINT Num1$; " = ";
      COLOR 2: PRINT USING "\    \"; Trim$((b));

      COLOR 15: LOCATE , 41: PRINT "A";
      COLOR 7: PRINT Num1$; " = ";
      COLOR 3: PRINT USING "\    \"; Trim$(SQR(a * a + b * b));
      COLOR 15: LOCATE , 58: PRINT "Phi";
      COLOR 7: PRINT Num1$; " = ";
      COLOR 4: PRINT USING "\    \"; Trim$(-ArcTan((b), (a)));
   NEXT n

   WHILE INKEY$ = "": WEND

END SUB

SUB ShowPage (PageNum AS INTEGER)

   ' Switches the visual page to PageNum, puts the active page in the
   ' background and clears it.

   SCREEN ScreenMode, 0, 1 - PageNum, PageNum
   LINE (0, 0)-(XSize - 1, YSize - 1), 0, BF

END SUB

' The Sinc function.
' Sinc(x) is equal to SIN(Pi * x) / (Pi * x).
FUNCTION Sinc# (x AS DOUBLE)

   IF x = 0 THEN
      Sinc = 1
   ELSE
      Sinc = SIN(Pi * x) / (Pi * x)
   END IF

END FUNCTION

' A building block function to make more complicated functions out of.
' The wave has period x, and duty cycle of DC.
' (Duty cycle = percentage of time the function is 1.)
FUNCTION SquareWave# (x AS DOUBLE, DC AS DOUBLE)

   NewX = Mod2(x, 1)

   IF NewX > 1 / 2 THEN NewX = NewX - 1

   IF (DC * 1 / 2) - ABS(NewX) >= 0 THEN
      SquareWave = 1
   ELSE
      SquareWave = 0
   END IF

END FUNCTION

' This function is used 4 times by GraphFourierCoeffs() to draw the 4 little
' graphs.
SUB SubGraph (d() AS DOUBLE, x0 AS INTEGER, y0 AS INTEGER, Label AS STRING)

' We have a screen space of 250x200 pixels, with top left corner (x0, y0).

   COLOR 9

   YMin = 30000: YMax = -30000
   FOR i = 0 TO 50
      IF d(i) < YMin THEN YMin = d(i)
      IF d(i) > YMax THEN YMax = d(i)
   NEXT i
   IF YMin > 0 THEN YMin = 0 ELSE IF YMax < 0 THEN YMax = 0
   YRange = YMax - YMin
   RYMin = YMin: RYMax = YMax: RYRange = YRange
   YMax = YMax + YRange * ExtraYSpace
   YMin = YMin - YRange * ExtraYSpace
   YRange = YMax - YMin
   IF YRange < .001 THEN YRange = 0

   ' Locate where in screen space we have y = 0, draw x axis.
   IF YRange = 0 THEN
      BaseY = 100
   ELSE
      BaseY = 200 + YMin / YRange * 200
   END IF
   IF BaseY >= 0 AND BaseY <= 200 THEN
      LINE (x0, y0 + BaseY)-(x0 + 252, y0 + BaseY), 1
   ELSE
      STOP
   END IF

   ' draw y axis and label
   LINE (x0, y0)-(x0, y0 + 200), 1
   LOCATE (y0 - 8) / 16, (x0 - 8) / 8: PRINT Label$;
   ' draw yMax, yMin and show on axis
   IF YRange = 0 THEN
      LINE (x0 - 2, y0 + BaseY)-(x0 + 2, y0 + BaseY), 1
      LOCATE (8 + y0 + BaseY) / 16, (x0 - 45) / 8
      PRINT USING "\   \"; Trim$(0);
   ELSE
      LINE (x0 - 2, y0 + BaseY - RYMin / YRange * 200)-(x0 + 2, y0 + BaseY - RYMin / YRange * 200), 1
      LOCATE (8 + y0 + BaseY - RYMin / YRange * 200) / 16, (x0 - 45) / 8
      PRINT USING "\   \"; Trim$((RYMin));
      LINE (x0 - 2, y0 + BaseY - RYMax / YRange * 200)-(x0 + 2, y0 + BaseY - RYMax / YRange * 200), 1
      LOCATE (8 + y0 + BaseY - RYMax / YRange * 200) / 16, (x0 - 45) / 8
      PRINT USING "\   \"; Trim$((RYMax));
   END IF

   ' graph data
   FOR n = 0 TO 50
      IF YRange = 0 THEN
         CIRCLE (x0 + n * 5, y0 + BaseY), 1, 14
      ELSE
         LINE (x0 + n * 5, y0 + BaseY)-(x0 + n * 5, y0 + BaseY - d(n) / YRange * 200), 14
         CIRCLE (x0 + n * 5, y0 + BaseY - d(n) / YRange * 200), 1, 14
      END IF
   NEXT n

END SUB

' Rounds of a number so it can fit nicely in a table of coefficients,
' and converts it into a string at the same time.
FUNCTION Trim$ (a AS DOUBLE)

   DIM x AS SINGLE
  
   x = INT(a * 10000) / 10000
   IF ABS(x) < .0001 THEN
      x = 0
   END IF

   Trim = LTRIM$(RTRIM$(STR$(x)))
  
END FUNCTION