// Fourier series routines #include #define MAXITER 15 double om; double PI = 3.1415926535; int flag; // analytical function to study: // ---------------------------- // Note: Exact Fourier coefficients for this periodic function are: // an=0 if n is even, -1/2n if n=4p+1, 1/2n if n=4p+3 // bn=0 if n=4p, -1/2n if n=4p+1, 1/n if n=4p+2, -1/2n if n=4p+3 //---------------------------------------------------------------- double F(double x) { if (x<-PI) return 0.0; else if (x<-PI/2) return PI/4; else if (x<=PI) return -PI/4; else return 0.0; } // Function to integrate by Romberg method double FUNC(double x) { if (!flag) return (F(x)*cos(om*x)); //for an else return (F(x)*sin(om*x)); //for bn } /*************************************************************** * Calculate the Fourier harmonic #n of a periodic discreet * * function F(x) defined by ndata points. * * ------------------------------------------------------------ * * Inputs: * * ndata: 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]. * * * * C++ 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 ----------------------------------------------------------------*/ void DiscreetFourierHn(int ndata, double *X, double *Y, int n, double *a, double *b) { int i; double xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp; T=X[ndata]-X[1]; xi=(X[ndata]+X[1])/2; om = 2*PI*n/T; *a=0; *b=0; for (i=1; iMAXITER) itermax=MAXITER; r = FUNC(a); ta = (r + FUNC(b)) / 2; *n=0; pas=b-a; t[0][0]=ta*pas; do { *n = *n + 1; pas = pas/2; s=ta; for (i=1; i=itermax) return t[*n][0]; } while (*obtprec >prec || *n