Representation of Fourier series depends on tabulation points

Question

Well, I had task to create function that does Fourier series with some mathematical function, so I found all the formulas, but the main problem is when I change count of point on some interval to draw those series I have very strange artifact:

This is Fourier series of sin(x) on interavl (-3.14; 314) with 100 point for tabulation

And this is same function with same interval but with 100000 points for tabulation

Code for Fourier series coeficients:

void fourieSeriesDecompose(std::function<double(double)> func, double period, long int iterations, double *&aParams, double *&bParams){
  aParams = new double[iterations];
  aParams[0] = integrateRiemans(func, 0, period, 1000);

  for(int i = 1; i < iterations; i++){
    auto sineFunc = [&](double x) -> double { return 2 * (func(x) * cos((2 * x * i * M_PI) / period)); };
    aParams[i] = integrateRiemans(sineFunc, -period / 2, period / 2, 1000) / period;
  }

  bParams = new double[iterations];
  for(int i = 1; i < iterations; i++){
      auto sineFunc = [&](double x) -> double { return  2 * (func(x) * sin(2 * (x * (i + 1) * M_PI) / period)); };
      bParams[i] = integrateRiemans(sineFunc, -period / 2, period / 2, 1000) / period;
  }

}

This code I use to reproduce function using found coeficients:

double fourieSeriesCompose(double x, double period, long iterations, double *aParams, double *bParams){
  double y = aParams[0];

  for(int i = 1; i < iterations; i++){
    y += sqrt(aParams[i] * aParams[i] + bParams[i] * bParams[i]) * cos((2 * i * x * M_PI) / period - atan(bParams[i] / aParams[i]));
  }
  return y;
}

And the runner code

double period = M_PI * 2;
auto startFunc = [](double x) -> double{ return sin(x); };

fourieSeriesDecompose(*startFunc, period, 1000, aCoeficients, bCoeficients);
auto readyFunc = [&](double x) -> double{ return fourieSeriesCompose(x, period, 1000, aCoeficients, bCoeficients); };

tabulateFunc(readyFunc);
scaleFunc();
//Draw methods after this

In all those loops where you multiply by `i` (the iteration number), you never divide by `iterations`, so larger numbers of iterations create larger and larger results, not more refined results. Was this the intent? — 1201ProgramAlarm, Dec 16 '20 at 23:00
@1201ProgramAlarm, I'll try to do this with division, but I haven't seen such thing in formulas — Psyhich, Dec 16 '20 at 23:05
@Psyhich: I think 1201 has a point. Think back about the basics of integration - you obtain the area under the curve by breaking it up in small rectangular strips of width `dx`, with `dx->0` in theory. Here, you set `dx` to `period/iterations`. So the area of each rectangular strip is `func(x) * dx` = `func(x) * period/iterations`. You need to divide. — MSalters, Dec 17 '20 at 02:29

Spektre · Answer 1 · 2020-12-17T09:23:29.470

see:

How to compute Discrete Fourier Transform?

So if I deciphered it correctly the aParams,bParams represent the real and imaginary part of the result then the angles in sin and cos must be the same but you have different! You got this:

auto sineFunc = [&](double x) -> double { return 2*(func(x)*cos((2* x* i   *M_PI)/period));
auto sineFunc = [&](double x) -> double { return 2*(func(x)*sin( 2*(x*(i+1)*M_PI)/period));

as you can see its not the same angle. Also what is period? You got iterations! if it is period of the function you want to transform then it should be applied to it and not to the kernel ... Also integrateRiemans does what? its the nested for loop to integrate the furrier transform? Btw. hope that func is real domain otherwise the integration/sumation needs both real and imaginary part not just one ...

So what you should do is:

create (cplx) table of the func(x) data on the interval you want with iterations samples

so for loop where x = x0+i*(x1-x0)/(iterations-1) and x0,x1 is the range you want the func to sample. Lets call it f[i]
```
for (i=0;i<iteration;i++) f[i]=func(x0+i*(x1-x0)/(iterations-1));
```
furrier transform it

something like this:
```
for (i=0;i<iteration;i++) a[i]=b[i]=0;
for (j=0;j<iteration;j++)
 for (i=0;i<iteration;i++)
  {
  a[j]+=f[i]*cos(-2.0*M_PI*i*j/iterations);
  b[j]+=f[i]*sin(-2.0*M_PI*i*j/iterations);
  }
```
now a[],b[] should hold your slow DFT result ... beware integer rounding ... depending on compiler you might need to cast some stuff to double to avoid integer rounding.

Representation of Fourier series depends on tabulation points

1 Answers1