C++

Forum

 
 
Fourier Series function plotting differently than graph calc

SophiaCristina (16)
I'm new on this forum and kinda new in c++ and certain maths...
I'm not sure, but I think I did a correct Fourier formula on Geogebra graph calculator, but, when using the same logic on c++, this is not working...
Graphs in Geogebra and my plot: https://i.imgur.com/pETifd4.png

But maybe my logic is wrong, but I still can't find where. I'm almost sure its a order logic mistake...
Sorry if this don't belongs here or another thing, english is not my native language and i though this was the most appropriate place to create this topic. The "Preview" button of this forum didn't work for me...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
 
  // FOURIER SERIES:
	void FourierSeries()
	{
		// Integra:
		int n2 = 30, TR, TG, TB;
		double a0 = 2 * IntgrlRAM(0, 1, 0, 1000, 1);
		double DDiv = Tau / 1000;
		cout << "a0 = " << a0 << endl;

		for (float tx = 0; tx <= Tau; tx = tx + 0.01)
		{
			double FS = 0, Sum = 0;
			for (int n = 0; n <= n2; ++n)
			{
				double LRAM = 0, RRAM = 0, a = 0, LRAM2 = 0, RRAM2 = 0, b = 0, aR = 0, bR = 0, aL = 0, bL = 0, as = 0, bs = 0;
				for (int m = 1; m <= 1000; ++m)
				{
					// FORMULA:
					// INTEGRAIS:
					aR = (sin(m * DDiv) * cos(Pi * n * (m * DDiv))) * DDiv;
					bR = (sin(m * DDiv) * sin(Pi * n * (m * DDiv))) * DDiv;
					aL = (sin(m * DDiv) * cos(Pi * (n - 1) * (m * DDiv))) * DDiv;
					bL = (sin(m * DDiv) * sin(Pi * (n - 1) * (m * DDiv))) * DDiv;
					RRAM = RRAM + aR;
					LRAM = LRAM + aL;
					a = ((RRAM - LRAM) * 0.5) + LRAM;
					RRAM2 = RRAM2 + bR;
					LRAM2 = LRAM2 + bL;
					b = ((RRAM2 - LRAM2) * 0.5) + LRAM2;
					if (n == 0) { a = a0; b = 0; }
					// #######
				}
				as = 2 * a * cos(Tau * n * tx);
				bs = 2 * b * sin(Tau * n * tx);
				Sum = Sum + (as + bs);
				cout << " | Sum:" << Sum;
			}

			FS = 0.5 * a0 + Sum;
			int xg = round((tx / Tau) * GradeX - 1) + BordaY + 1;
			int yg = (GradeY - round((FS / 2.58759) * (GradeY - 2))) + BordaX + 1;
			if (InImg(xg, yg))
			{
			LinearRGB((tx / Tau), 1, 1, TR, TG, TB, false);
			unsigned char Clr[] = { TR, TG, TB };
			Imagem.draw_point(xg, yg, Clr);
			}
		}
	}


Notes:
* "LinearRGB" is just a function that by reference gives the RGB based on a number between 0.0 and 1.0, so ignore;
* xg and yg are just variables for the plotter, ignore;
* "Clr[]" and "Imagem" are CImg classes, ignore;
* The cout prints a lot of lines, you may take it off;
* The division number (2.58759) on "yg" is based on the number i got from this lot of cout lines.

"IntgrlRAM(0, 1, 0, 1000, 1)" is my Integral Retangle Aproximation Method, the code is this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
// INTEGRAL:
double IntgrlRAM(double a, double b, int m, int Iteracao, short int OpRet)
{
	if (OpRet > 3) { OpRet = 3; } if (OpRet < 1) { OpRet = 1; }
	if (a > b) { double Tmp; Tmp = a; a = b; b = Tmp; }
	double LRAM = 0, RRAM = 0, MRAM = 0;
	double DDiv = (b - a) / Iteracao;

	for (int n = 1; n <= Iteracao; ++n)
	{
		// FORMULA:
		RRAM = RRAM + (MiniForm((n * DDiv) + a, 1) * DDiv);
		LRAM = LRAM + (MiniForm(((n - 1) * DDiv) + a, 1) * DDiv);
		MRAM = ((RRAM - LRAM) * 0.5) + LRAM;
		// #######
	}
	if (OpRet == 1) { return(MRAM); } if (OpRet == 2) { return(LRAM); } if (OpRet == 3) { return(RRAM); }
}


* m is just a variable i use when i change the formula...
* "MiniForm(n, Frequency)" is just a f(x) form, its set on "sin(n)" in this case...
jonnin (11333)
just making sure, are you aware that c++ trig functions are in radians and most calculators default to degrees (but can be set to rads)? Is that a possible issue?
lastchance (6980)
What function are you trying to find a Fourier series for?

Your loop from lines 16 to 32 is very obscure, but it looks as if it is trying to do two very-long-winded integrations by trapezium rule for the sine and cosine parts of the series. But surely IntgrlRAM() is already set up to do this, so why aren't you using it?

Also your a and b coefficients for a Fourier series are fixed constants. You should work them out just once. Yet your code appears to be working them out repeatedly for successive tx values.
Last edited on
SophiaCristina (16)
@jonnin
Thanks, but Geogebra, when you put a number directly, its degrees, but when the number is a function of something, its radians, basically, when you put a number directly, it appends a "º" after your number, but, deleting this "º" it turns to radians, so, if you write sin(x) or sin(n), its radians, but if you write sin(3.14), it morphs to sin(3.14º), but if you delete "º", it returns the radians value.
I investigated it anyway, just to be sure, but looks fine after all investigation... :/

@lastchance
Resuming, i do musics and i'm learning maths and i'm now into this Fourier things, its useful for me, since I can make square waves from sines (but i already got a summation formula for it), I went to wikipedia and used the first formula that appears:
https://en.wikipedia.org/wiki/Fourier_series#Definition
Where the formula is, its until before section "A more general definition"...
The function s(x) is just a sin(x) on my program, that is the one i'm trying to find a fourier for, but its for learning, i'll later use another formulas...

The 16 - 32, i made because, i'm noob on c++ and i have no parsing lib yet, so, the IntgrlRAM() is integrating the "sin(x)" only, but 16 - 32 is integrating "sin(x) * cos(2 * pi * n * x)", to the integral have no conflicts with "tx" in its approximation, i made the "x" in this formula the "(m * DDiv)".

Hmm, i declared "a" and "b" in all different brackets and without "= 0" and the result is the same... :/
Sorry, i understood your reply of "fixed const" as the idea of declaring it as "= 0", sorry if i interpreted it wrong and about my "layness".
I'm noob in maths, and my idea was to use this for-tx to bring a value for each x in the formula "as = 2 * a * cos(Tau * n * tx); bs = 2 * b * sin(Tau * n * tx);".
I'll try to find a value for a single x, then i reply this topic when i do it...

EDIT: Derp, why trying to find the value of a single x?! Its the same as tx, but returning a single value, anyway, this is the print of the cout in "x = Pi":
https://imgur.com/YHafBQf

---

Anyway, thanks for the replies...
Last edited on
lastchance (6980)
You can have a play around with Fourier series in the following code (just write a function to use in the constructor called from main()).

An example for a square-wave approximation using 20 harmonics:
https://imgur.com/a/RHMWVkU

You can just redirect the output to a file and plot with gnuplot.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
 
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <functional>
#include <algorithm>
#include <cmath>
using namespace std;

const double PI = 3.14159265358979;

using func = function<double(double)>;                     // Function object

double integral( func f, double a, double b, int n );


//==========================================================


class FourierSeries
{
   double period;                                          // Period
   double omega;                                           // Fundamental frequency
   int nmax;                                               // Highest harmonic
   vector<double> a, b;
public:
   FourierSeries( double (*f)( double ), double T, int N ) : period( T ), nmax( N )
   {
      omega = 2.0 * PI / period;          
      a = vector<double>( 1 + nmax, 0 );   b = a;
      int Ndt = max( 10 * N, 100 );                        // Number of intervals for numerical integration
      double factor = 2.0 / period;                        // Normalising factor for basis functions

      for ( int n = 0; n <= nmax; n++ )                    // Compute coefficients
      {
         func c = [=]( double t ){ return f( t ) * cos( n * omega * t ); };
         func s = [=]( double t ){ return f( t ) * sin( n * omega * t ); };
         a[n] = factor * integral( c, 0.0, period, Ndt );
         b[n] = factor * integral( s, 0.0, period, Ndt );
      }
   }

   double evaluate( double t )                                       // Sum Fourier series
   {
      double sum = 0.5 * a[0];
      for ( int n = 1; n <= nmax; n++ ) sum += a[n] * cos( n * omega * t ) + b[n] * sin( n * omega * t );
      return sum;
   }

   void tabulate( ostream &out, double t0, double dt, int nt )       // Output
   {
      #define FMT << fixed << setprecision( 4 ) << setw( 15 ) <<
      for ( int i = 0; i <= nt; i++ )
      {
         double t = t0 + i * dt;
         out FMT t FMT evaluate( t ) << '\n';
      }
   }
};


//==========================================================


double integral( func f, double a, double b, int n )                 // Numerical integration: mid-ordinate rule
{
   double dt = ( b - a ) / n;
   double sum = 0;
   for ( int i = 1; i <= n; i++ ) sum += f( a + ( i - 0.5 ) * dt) * dt;
   return sum;
}


//==========================================================


double squareWave( double t )
{  
   t = t - floor( t );
   return t < 0.5 ? 1 : 0;
}


double sawTooth( double t )
{
   t = t - floor( t );
   return t;
}


//==========================================================


int main()
{
   const double period = 1.0;                    // Period of function
   const int N = 20;                             // Highest harmonic

   FourierSeries fs( squareWave, period, N );    // Try other functions besides squareWave
   fs.tabulate( cout, 0.0, period / 100, 400 );
}


//==========================================================
Last edited on
SophiaCristina (16)
Oh, you made the code for me... :)
Thank you very much for the effort... Now i can see i have a lot to learn in c++ yet...

Here is the plotted result with CImg in a bad scale:
https://imgur.com/almqM66

By the way, I don't knew about gnuplot, for math i was using a software called SageMath, i'll research about it.
Last edited on
Topic archived. No new replies allowed.