From db2d4333efc5af96d8cddea245463198208d6043 Mon Sep 17 00:00:00 2001
From: koldo <koldo@f0d560ea-af0d-0410-9eb7-867de7ffcac7>
Date: Fri, 12 Feb 2016 22:24:24 +0000
Subject: [PATCH] Eigen_demo: Added FFT demo

git-svn-id: svn://ultimatepp.org/upp/trunk@9470 f0d560ea-af0d-0410-9eb7-867de7ffcac7
---
 reference/Eigen_demo/Eigen_demo.upp |  3 +-
 reference/Eigen_demo/eigen_demo.cpp | 10 ++++
 reference/Eigen_demo/fft.cpp        | 84 +++++++++++++++++++++++++++++
 3 files changed, 96 insertions(+), 1 deletion(-)
 create mode 100644 reference/Eigen_demo/fft.cpp

diff --git a/reference/Eigen_demo/Eigen_demo.upp b/reference/Eigen_demo/Eigen_demo.upp
index 1eb09f6cd..d39ae16a9 100644
--- a/reference/Eigen_demo/Eigen_demo.upp
+++ b/reference/Eigen_demo/Eigen_demo.upp
@@ -6,7 +6,8 @@ uses
 
 file
 	eigen_demo.cpp,
-	non-linear.cpp;
+	non-linear.cpp,
+	fft.cpp;
 
 mainconfig
 	"" = "";
diff --git a/reference/Eigen_demo/eigen_demo.cpp b/reference/Eigen_demo/eigen_demo.cpp
index 03296252a..2874563d8 100644
--- a/reference/Eigen_demo/eigen_demo.cpp
+++ b/reference/Eigen_demo/eigen_demo.cpp
@@ -7,6 +7,7 @@ using namespace Upp;
 using namespace Eigen;
 
 void NonLinearTests();
+void FFTTests();
 
 struct SerialTest {
 	MatrixXd m;
@@ -30,6 +31,10 @@ struct SerialTest {
 		
 CONSOLE_APP_MAIN
 {	
+	FFTTests();
+	
+	Cout() << "\nFin";
+	getchar();
 	// http://eigen.tuxfamily.org/dox/TutorialMatrixClass.html
 	Cout() << "\n\nTutorial page 1 - The Matrix class";
 	
@@ -695,6 +700,11 @@ CONSOLE_APP_MAIN
 		
 	NonLinearTests();
 	
+	Cout() << "\nPress enter to continue\n";
+	ReadStdIn();
+	
+	FFTTests();
+	
 	Cout() << "\nPress enter to end";
 	ReadStdIn();
 }
diff --git a/reference/Eigen_demo/fft.cpp b/reference/Eigen_demo/fft.cpp
new file mode 100644
index 000000000..dd72c068e
--- /dev/null
+++ b/reference/Eigen_demo/fft.cpp
@@ -0,0 +1,84 @@
+#include <Core/Core.h>
+
+using namespace Upp;
+
+#include <plugin/Eigen/Eigen.h>
+
+#ifdef USE_FFTW
+#include <fftw3.h>
+#endif
+
+using namespace Eigen;
+
+
+void FFTTests()
+{
+	Cout() << "\n\nFFT sample\nGets the FFT of equation"
+			  "\n f(t) = 2*sin(2*PI*t/50 - PI/3) + 5*sin(2*PI*t/30 - PI/2) + 30*sin(2*PI*t/10 - PI/5)"
+			  "\nsampled with a frequency of 14 samples/second";
+	
+	int numData = 8000;
+	double samplingFrecuency = 14;	
+	String csvSep = ";";
+
+    // Filling the data series
+    VectorXd timebuf(numData);
+    {
+	    double t = 0;
+	    for (int i = 0; i < numData; ++i, t = i*1/samplingFrecuency) 
+	       	timebuf[i] = 2*sin(2*M_PI*t/50 - M_PI/3) + 5*sin(2*M_PI*t/30 - M_PI/2) + 30*sin(2*M_PI*t/10 - M_PI/5);
+    }
+    
+    // FFT
+    VectorXcd freqbuf;
+    FFT<double> fft;
+    fft.SetFlag(fft.HalfSpectrum);
+    fft.fwd(freqbuf, timebuf);
+	
+	// Filter the FFT. Frequemcies between 1/25 and 1/35 Hz are removed
+	// Original FFT is not changed for saving it later
+	VectorXcd freqbuf2(freqbuf.size());
+	{
+	    for (int i = 0; i < freqbuf.size(); ++i) {
+	        double freq = i*samplingFrecuency/numData;
+	        double T = 1/freq;
+	        if (T > 25 && T < 35)
+	            freqbuf2[i] = 0;
+	        else
+	            freqbuf2[i] = freqbuf[i];
+	    }
+	}
+	
+	// Inverse filtered FFT to get filtered series
+	VectorXd timebuf2(numData);
+	fft.inv(timebuf2, freqbuf2);
+	
+	// Saving original and filtered FFT
+	{
+	    String str;
+	    str << "Frec" << csvSep << "T" << csvSep << "fft" << csvSep << "Filtered fft";
+	    for (int i = 0; i < freqbuf.size(); ++i) {
+	        double freq = i*samplingFrecuency/numData;
+	        double T = 1/freq;
+	        str << "\n" << freq << csvSep << (freq > 0 ? FormatDouble(T) : "") << csvSep 
+	        			<< 2*std::abs(freqbuf[i])/numData << csvSep 
+	        			<< 2*std::abs(freqbuf2[i])/numData;;
+	    }
+	    String fftFileName = AppendFileName(GetDesktopFolder(), "fft.csv");
+	    Cout() << "\nFFT saved in '" << fftFileName << "'";
+	    SaveFile(fftFileName, str);
+	}
+	
+	// Saving original and filtered series
+	{
+	    String str;
+	    str << "Time" << csvSep << "Data" << csvSep << "Filtered data";
+	    double t = 0;
+	    for (int i = 0; i < numData; ++i, t = i*1/samplingFrecuency) 
+	       	str << "\n" << t << csvSep << timebuf[i] << csvSep << timebuf2[i];;
+	    String dataFileName = AppendFileName(GetDesktopFolder(), "data.csv");
+	    Cout() << "\nSource data saved in '" << dataFileName << "'";
+	    SaveFile(dataFileName, str);
+    }
+}
+