fourier, readme

1 files changed, 82 insertions, 0 deletions

diff --git a/discrete_fourier.c b/discrete_fourier.c
new file mode 100644
index 0000000..68a5d98
--- /dev/null
+++ b/discrete_fourier.c
@@ -0,0 +1,82 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#define _PI 3.14159
+
+#define _N 100
+
+struct ret
+{
+	double *val;
+	unsigned int valc;
+};
+
+static double absolute(double a, double b);
+
+struct ret dft(unsigned long N, double *val);
+
+int main(void)
+{
+  double val[_N];
+  for(int i = 0; i < _N ; i++)
+  {
+    val[i] = 0;
+  }
+
+  for(int i = 0; i < _N; i+=4)
+  {
+    val[i] = 1;
+  }
+
+  struct ret a = dft(_N, (double*)&val);
+
+  for(int i = 0; i < a.valc; i++)
+  {
+    printf("%f,",a.val[i]);
+  }
+
+  printf("\n");
+
+  return 0;
+}
+
+//discrete fourier transform
+// N: Number of values
+// *val: Array of values
+struct ret dft(unsigned long N, double *val)
+{
+	//Dynamic memory allocation
+	struct ret _ret;
+	double kmax = N / 2;
+	double *ival = malloc(sizeof(double) * kmax); //Array of imaginary values
+	double *rval = malloc(sizeof(double) * kmax); //Array of real values
+
+	_ret.val = malloc(sizeof(double) * kmax);
+	_ret.valc = kmax;
+
+	//Calculate Fourier-Transform for every k
+	for (unsigned int k = 0; k < kmax; k++)
+	{
+		rval[k] = 0;
+		ival[k] = 0;
+		//Actual discrete Fourier-Transform
+		for (unsigned int n = 0; n < N; n++)
+		{
+			double angle = (2 * _PI * k * n) / N;
+			double rx =   val[n] * cos(angle);//real part
+			double ix = - val[n] * sin(angle);//imaginary part
+
+			rval[k] += rx;
+			ival[k] += ix;
+		}
+		_ret.val[k] = absolute(rval[k], ival[k]);
+	}
+
+	return _ret;
+}
+
+static double absolute(double a, double b)
+{
+	return sqrt(pow(a, 2) + pow(b, 2));
+}