'********************************************************** '* 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