- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

F90 Version By J-P Moreau, Paris. (www.jpmoreau.fr) published a method to calculate the coefficients a and b for the standard Fourier series.

I am taking the a and b symbols from the Kaplan and Lewis textbook.

The code is

```
! Fourier series routines (version with dynamic allocations)
! analytical function(s) 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
! ----------------------------------------------------------------
REAL*8 Function F(x) !example 1
real*8 x, PI
PI=4.d0*datan(1.d0)
if (x<-PI) then
F = 0.d0
else if (x<-PI/2.d0) then
F = PI/4.d0
else if (x<=PI) then
F = -PI/4.d0
else
F = 0.d0
end if
END
! Note: Exact Fourier coefficients for this periodic function are:
! a0=1/pi, an = 0 if n is odd, an = -2/(pi*(n2-1)) if n is even;
! b1=1/2, bn = 0 for n>=2.
!-----------------------------------------------------------------
REAL*8 Function F1(x) !example 2
real*8 x, PI
PI=4.d0*datan(1.d0)
if (x<0.d0) then
F1=0.d0
else if (x<PI) then
F1=dsin(x)
else
F1=0.d0
end if
END
! Function to integrate by Romberg method
REAL*8 Function FUNC(om,x,flag)
real*8 om,x; integer flag
real*8 F, F1
if (flag==0) then
FUNC = (F1(x)*dcos(om*x)) !for an
else
FUNC = (F1(x)*dsin(om*x)) !for bn
end if
END
!****************************************************************
!* 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]. *
!* *
!* F90 Version By J-P Moreau, Paris. *
!* (www.jpmoreau.fr) *
!****************************************************************
! 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(ndata, X, Y, n, a, b)
real*8 X(1:ndata), Y(1:ndata), a, b
integer ndata,n,i
real*8 xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp
real*8 PI
PI=4.d0*datan(1.d0)
T=X(ndata)-X(1); $xd[$ndata] - $xd[1]
om = 2*PI*n/T; a=0.d0; b=0.d0
do i=1, ndata-1
wa=X(i); wb=X(i+1)
wc=Y(i); wd=Y(i+1)
if (wa.ne.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.d0*(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
end do
a = a/T; b = b/T
if (n.ne.0) then
a = a*2.d0/om; b = b*2.d0/om
end if
return
End
!*******************************************************
!* Integral of a function FUNC(X) by Romberg's method *
!* --------------------------------------------------- *
!* INPUTS: *
!* a begin value of x variable *
!* b end value of x variable *
!* prec desired precision *
!* itermin minimum number of iterations *
!* itermax maximum number of iterations *
!* *
!* OUTPUTS: *
!* obtprec obtained precision for integral *
!* niter number of iterations done *
!* integral the integral of FUNC(X) from a to b *
!* *
!*******************************************************
real*8 Function RombergIntegral(a,b,om,flag,prec,obtprec,niter,itermin,itermax)
parameter(MAXITER = 15)
integer flag
real*8 a,b,om,prec,obtprec,FUNC
real*8 t(MAXITER,MAXITER)
real*8 pas,r,s,ta
if (itermax>MAXITER) itermax=MAXITER
r = FUNC(om,a,flag)
ta = (r + FUNC(om,b,flag)) / 2.d0
niter=0
pas=b-a
t(0,0)=ta*pas
100 niter=niter+1
pas=pas/2.d0
s=ta
do i=1, 2**niter-1
s = s + FUNC(om,a+pas*i,flag)
end do
t(0,niter)=s*pas
r=1.d0
do i=1, niter
r=r*4.d0
j=niter-i
t(i,j)=(r*t(i-1,j+1) - t(i-1,j))/(r-1.d0)
end do
obtprec = DABS(t(niter,0) - t(niter-1,0))
if (niter > itermax) goto 200
if (niter<itermin) goto 100
if (obtprec>prec) goto 100
200 RombergIntegral = t(niter,0)
return
end
!****************************************************************
!* Calculate the Fourier harmonic #n of a periodic function F(x)*
!* analytically defined. *
!* ------------------------------------------------------------ *
!* Inputs: *
!* t1 : -period/2 *
!* t2 : period/2 *
!* n : order of harmonic *
!* *
!* 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]. *
!* *
!* F90 Version By J-P Moreau, Paris. *
!* (www.jpmoreau.fr) *
!****************************************************************
! Note: When a periodic function f(x), of period T, can be developped
! in a Fourier series, this one is unique:
! n=inf.
! F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)
! n=1
! (b0 is always zero).!
!
! The coefficients an, bn are given by:
! T/2
! a0 = 1/T Sum (f(x) dx)
! -T/2
! T/2
! an = 2/T Sum (f(x) cos(2npi/T x) dx)
! -T/2
! T/2
! an = 2/T Sum (f(x) sin(2npi/T x) dx)
! -T/2
! Here, the integrals are calculated by the Romberg method.
! --------------------------------------------------------------------
Subroutine AnalyticFourierHn(t1,t2,n,a,b)
real*8 t1,t2,a,b,om
real*8 T,precision,temp, PI
integer flag
PI=4.d0*datan(1.d0)
T=dabs(t2-t1)
om=2.d0*PI*n/T
precision=1.d-10; itmax=15
flag=0 ! calculate an
a=RombergIntegral(t1,t2,om,flag,precision,temp,niter,5,itmax)*2.d0/T
flag=1 ! calculate bn
b=RombergIntegral(t1,t2,om,flag,precision,temp,niter,5,itmax)*2.d0/T
if (n==0) a=a/2.d0
return
End
! end of file fourier.f90
```

part of a set of functions. I have some data in a mySQL database in the cloud and I was wanting to do a Fourier function for the temperature, for some statistical analysis.

So it has been an interesting challenge to translate the Fortran into PHP. PHP is not great with errors, so you go one line at a time, Dreamweaver is not great with formatting automatically , so the do and if statements over a long set of code can become misaligned.

Anyway I got the code running in PHP but matching the input data

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

was proving to be challenging, PHP is not nice to debug with the long cycle times and the lack of nice error output, one ends up using a lot of log files.

So to check the code I created a SLN file for the Fortran and used the same data from the php. The sln file is enclosed.

The final PHP output

and the Fortran output are the same - as expected,

But the Fourier equation has an inversion in the summation, I need to use - instead of +.

```
do j = 1,i
z(j) = a(1) - ((a(2)*cos((2.0*3.142/float(i-1))*j)) + (b(2)*sin((2.0*3.142/float(i-1))*j)))
write(14,200)j,x(j),y(j),z(j)
200 Format(i5,3(" , ",f10.3))
```

The same is true for the PHP.

I will now sit down with my numerical analysis textbook and work out what the problem is in the code, but the French code has the same error in all languages.

Any way it is interesting.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

In the comments for subroutine DiscreetFourierHn it lists a and b are outputs, however the code performs summations into a and b with a and b on rhs. IOW the code assumes a and b are used as inputs as well. Either the comment is wrong or the caller must initialize (?0.0, n.m?) these values prior to call.

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

a and b are initialized to zero in the original Fortran code and the PHP code before the summation starts.

The original book is from 1991 and is French, I will have to see if I can find it.

The Numerical Ideas for using this technique appear in the late 1950s, the best one is in the Acoustical Journal, but I need to go to the library to find it. You find the odd reference to the method, but usually it is PPT notes for a lecture where the Professor is demonstrating their math skills and ability to make a ham and cheese sandwich into some French dish.

I have started to back determine the algorithm from the code.

The first step is to transform the 1 -> N from -pi/2 to pi/2.

The next step is to determine the slope on the graph from i to i+ 1 using the standard Newtonian Formula,

I am looking at rewriting the a and b summations with standard symbols so they make some sense, the w's are hard to follow.

Plus I am only using n = 1 and so as long as the temp is a smooth curve it is ok, longer ones will need higher n's.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Kaplan and Lewis, Calculus and Linear Algebra - provides a sample calculation for the absolute function from -pi to +p for the Fourier series calculation.

The results from the integration are shown on the picture

The Fourier program predicts the a0, a1, etc to a reasonable number of decimal places, the a(1) value is very close, the first wave is shown

But the final equation still needs the form

`z(j) = a(1) - ((a(2)*cos((2.0*3.142/float(i-1))*j)) + (b(2)*sin((2.0*3.142/float(i-1))*j)))`

to get the correct answers, otherwise the answers are 1/10th of the real answer.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

@jimdempseyatthecove , it was my mistake.

`z(j) = a(1) + a(2)*cos(x(j)) + b(2)*(sin(x(j)))`

I got the Fourier equation wrong, the call is incorrect as you need to use the calculated T in the subroutine in the Fourier equation as he wrote the equation or use the X values directly if they are in Pi format.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page