'**********************************************************
'* This program calculates the Fourier coefficients of a  *
'* periodic discrete function F(x) by using the procedure *
'* DiscreetFourierHn of unit Fourier.                     *
'* ------------------------------------------------------ *
'* SAMPLE RUN:                                            *
'*                                                        *
'* F(x) has the following shape:                          *
'*                                                        *
'*  1 ---- 2  pi/4                                        *
'*    !  !                                                *
'*    !  !                                                *
'* ------!-------------------------------------------> x  *
'*  -pi  !    0         ! pi                              *
'*       !              !                                 *
'*     3 ---------------- 4 -pi/4                         *
'*                                                        *
'* Calculate the Fourier coefficients of a periodic       *
'* discrete function:                                     *
'*                                                        *
'* Number of points:  4                                   *
'*                                                        *
'*  x1 = -3.1415927                                       *
'*  y1 = 0.7853982                                        *
'*                                                        *
'*  x2 = -1.5707963                                       *
'*  y2 = 0.7853982                                        *
'*                                                        *
'*  x3 = -1.5707963                                       *
'*  y3 = -0.7853982                                       *
'*                                                        *
'*  x4 = 3.1415927                                        *
'*  y4 = -0.7853982                                       *
'*                                                        *
'* Lowest  harmonic: 0                                    *
'* Highest harmonic: 10                                   *
'*                                                        *
'* a0 = -0.392699087                                      *
'* b0 =  0.000000000                                      *
'*                                                        *
'* a1 = -0.500000023                                      *
'* b1 = -0.500000048                                      *
'*                                                        *
'* a2 = 0.0000000025                                      *
'* b2 = 0.5000000023                                      *
'* -----------------                                      *
'* a9 = -0.055555558                                      *
'* b9 = -0.055555583                                      *
'*                                                        *
'* a10 = -0.000000025                                     *
'* b10 =  0.100000005                                     *
'*                                                        *
'*                   Basic version by J-P Moreau, Paris.  *
'*                            (www.jpmoreau.fr)           *
'**********************************************************
DEFDBL A-H, O-Z
DEFINT I-N

'integer  icount,ih1,ih2,i,npoints
'double   a, b

DIM x(25), y(25)

  CLS
  PRINT
  PRINT " Calculate the Fourier coefficients of a periodic discrete function:"
  PRINT
  INPUT " Number of points: ", npoints
  PRINT

  FOR i = 1 TO npoints
    PRINT "  x"; i; " = "; : INPUT "", x(i)
    PRINT "  y"; i; " = "; : INPUT "", y(i)
	PRINT
  NEXT i
  INPUT " Lowest  harmonic: ", ih1
  INPUT " Highest harmonic: ", ih2
  PRINT
  icount = 1
  FOR i = ih1 TO ih2
    GOSUB 1000  'call DiscreetFourierHn(npoints,X,Y,i,a,b)
    PRINT "  a"; i; " = "; : PRINT USING "##.#########"; a
    PRINT "  b"; i; " = "; : PRINT USING "##.#########"; b
    IF (icount MOD 5) = 0 THEN
      icount = 0
      INPUT "", R$
    END IF
    icount = icount + 1
  NEXT i
  PRINT

END


'****************************************************************
'*   Calculate the Fourier harmonic #n of a periodic discreet   *
'*   function F(x) defined by npoints points.                   *
'* ------------------------------------------------------------ *
'* Inputs:                                                      *
'*            npoints: number of points of discreet function.   *
'*            X    : pointer to table storing xi abscissas.     *
'*            Y    : pointer to table storing yi ordinates.     *
'*                                                              *
'* Outputs:                                                     *
'*            a    : coefficient an of the Fourier series.      *
'*            b:   : coefficient bn of the Fourier series.      *
'* ------------------------------------------------------------ *
'* Reference: "Mathematiques en Turbo-Pascal By Marc Ducamp and *
'*             Alain Reverchon, 1. Analyse, Editions Eyrolles,  *
'*             Paris, 1991" [BIBLI 03].                         *
'*                                                              *
'*                         Basic Version By J-P Moreau, Paris.  *
'****************************************************************
' Note: The Fourier series  of a periodic discreet function F(x)
'       can be written under the form:
'                   n=inf.
'       F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)
'                   n=1
' ----------------------------------------------------------------
'Subroutine DiscreetFourierHn(npoints, X, Y, n, a, b)
1000 n = i
  pi = 3.1415926535#
  T = x(npoints) - x(1): xi = (x(npoints) + x(1)) / 2#
  om = 2# * pi * n / T: a = 0#: b = 0#
  FOR ii = 1 TO npoints - 1
    wa = x(ii): wb = x(ii + 1)
    wc = y(ii): wd = y(ii + 1)
    IF wa <> wb THEN
      wg = (wd - wc) / (wb - wa)
      wh = om * (wa - xi): wi = om * (wb - xi)
      IF n = 0 THEN
	a = a + (wb - wa) * (wc + wg / 2# * (wb - wa))
      ELSE
	wl = COS(wh): wm = SIN(wh)
	wn = COS(wi): wp = SIN(wi)
	a = a + wg / om * (wn - wl) + wd * wp - wc * wm
	b = b + wg / om * (wp - wm) - wd * wn + wc * wl
      END IF
    END IF
  NEXT ii
  a = a / T: b = b / T
  IF n <> 0 THEN
    a = a * 2# / om: b = b * 2# / om
  END IF
RETURN

' end of file discfour.bas