diff --git a/FDTD/excitation.cpp b/FDTD/excitation.cpp
index 290e254..5fc15e7 100644
--- a/FDTD/excitation.cpp
+++ b/FDTD/excitation.cpp
@@ -16,6 +16,7 @@
 */
 
 #include "tools/array_ops.h"
+#include "tools/useful.h"
 #include "fparser.hh"
 #include "tinyxml.h"
 #include "excitation.h"
@@ -100,13 +101,6 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS )
 	return true;
 }
 
-unsigned int Excitation::CalcNyquistNum(double fmax) const
-{
-	if (dT==0) return 1;
-	double T0 = 1/fmax;
-	return floor(T0/2/dT);
-}
-
 unsigned int Excitation::GetMaxExcitationTimestep() const
 {
 	FDTD_FLOAT maxAmp=0;
@@ -147,7 +141,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
 		Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0+fc) );
+	SetNyquistNum( CalcNyquistNum(f0+fc,dT) );
 }
 
 void Excitation::CalcDiracPulsExcitation()
@@ -216,7 +210,7 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
 		Signal_curr[n] = fParse.Eval(vars);
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0) );
+	SetNyquistNum( CalcNyquistNum(f0,dT) );
 }
 
 void Excitation::CalcSinusExcitation(double f0, int nTS)
@@ -240,7 +234,7 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
 		Signal_curr[n] = sin(2.0*PI*f0*t);
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0) );
+	SetNyquistNum( CalcNyquistNum(f0,dT) );
 }
 
 void Excitation::setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,
diff --git a/FDTD/excitation.h b/FDTD/excitation.h
index f35e5ba..bc31df2 100644
--- a/FDTD/excitation.h
+++ b/FDTD/excitation.h
@@ -34,7 +34,6 @@ public:
 
 	//! Get the excitation timestep with the (first) max amplitude
 	virtual unsigned int GetMaxExcitationTimestep() const;
-	unsigned int CalcNyquistNum(double fmax) const;
 
 	void setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,
 								 vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir );
diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp
index f0ea782..e0efb07 100644
--- a/FDTD/processcurrent.cpp
+++ b/FDTD/processcurrent.cpp
@@ -18,6 +18,7 @@
 #include "tools/global.h"
 #include "processcurrent.h"
 #include <iomanip>
+#include <complex.h>
 
 ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng)
 {
@@ -25,6 +26,7 @@ ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng)
 
 ProcessCurrent::~ProcessCurrent()
 {
+	ProcessCurrent::FlushData();
 }
 
 void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
@@ -44,6 +46,14 @@ void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
 	}
 }
 
+void ProcessCurrent::Init()
+{
+	FD_currents.clear();
+	for (size_t n=0;n<m_FD_Samples.size();++n)
+		FD_currents.push_back(0);
+}
+
+
 int ProcessCurrent::Process()
 {
 	if (Enabled==false) return -1;
@@ -142,11 +152,35 @@ int ProcessCurrent::Process()
 		break;
 	}
 
-//	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
+	//	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
 	current*=m_weight;
-	v_current.push_back(current);
-	//current is sampled half a timestep later then the voltages
-	file  << setprecision(m_precision) << (0.5 + (double)Eng->GetNumberOfTimesteps())*Op->GetTimestep() << "\t" << current << endl;
+
+	if (ProcessInterval)
+	{
+		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
+		{
+			v_current.push_back(current);
+			//current is sampled half a timestep later then the voltages
+			file  << setprecision(m_precision) << (0.5 + (double)Eng->GetNumberOfTimesteps())*Op->GetTimestep() << "\t" << current << endl;
+		}
+	}
+
+	if (m_FD_Interval)
+	{
+		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0)
+		{
+			double T = ((double)Eng->GetNumberOfTimesteps() + 0.5) * Op->GetTimestep();
+			for (size_t n=0;n<m_FD_Samples.size();++n)
+			{
+				FD_currents.at(n) += current * cexp( -2.0 * 1.0i * M_PI * m_FD_Samples.at(n) * T );
+			}
+			++m_FD_SampleCount;
+			if (m_Flush)
+					FlushData();
+			m_Flush = false;
+		}
+	}
+
 	return GetNextInterval();
 }
 
@@ -154,3 +188,8 @@ void ProcessCurrent::DumpBox2File( string vtkfilenameprefix, bool /*dualMesh*/ )
 {
 	Processing::DumpBox2File( vtkfilenameprefix, true );
 }
+
+void ProcessCurrent::FlushData()
+{
+	Dump_FD_Data(FD_currents,1.0/(double)m_FD_SampleCount,m_filename + "_FD");
+}
diff --git a/FDTD/processcurrent.h b/FDTD/processcurrent.h
index 247ec7d..571cf3c 100644
--- a/FDTD/processcurrent.h
+++ b/FDTD/processcurrent.h
@@ -28,11 +28,16 @@ public:
 
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
 
+	virtual void Init();
+	virtual void FlushData();
 	virtual int Process();
 
 	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
 protected:
 	vector<FDTD_FLOAT> v_current;
+
+	vector<_Complex double> FD_currents;
+	void WriteFDCurrents();
 };
 
 #endif // PROCESSCURRENT_H
diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp
index 931fcfa..1856670 100644
--- a/FDTD/processing.cpp
+++ b/FDTD/processing.cpp
@@ -16,7 +16,10 @@
 */
 
 #include "tools/global.h"
+#include "tools/useful.h"
 #include "processing.h"
+#include <complex.h>
+#include <climits>
 
 Processing::Processing(Operator* op, Engine* eng)
 {
@@ -26,7 +29,10 @@ Processing::Processing(Operator* op, Engine* eng)
 	m_PS_pos = 0;
 	SetPrecision(12);
 	ProcessInterval=0;
+	m_FD_SampleCount=0;
+	m_FD_Interval=0;
 	m_weight=1;
+	m_Flush = false;
 }
 
 Processing::~Processing()
@@ -53,21 +59,36 @@ bool Processing::CheckTimestep()
 	{
 		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) return true;
 	}
+
+	if (m_FD_Interval)
+	{
+		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0) return true;
+	}
 	return false;
 }
 
 int Processing::GetNextInterval() const
 {
 	if (Enabled==false) return -1;
-	unsigned int next=-1;
+	int next=INT_MAX;
 	if (m_ProcessSteps.size()>m_PS_pos)
 	{
-		next = m_ProcessSteps.at(m_PS_pos)-Eng->GetNumberOfTimesteps();
+		next = (int)m_ProcessSteps.at(m_PS_pos)-(int)Eng->GetNumberOfTimesteps();
+	}
+	if (ProcessInterval!=0)
+	{
+		int next_Interval = (int)ProcessInterval - (int)Eng->GetNumberOfTimesteps()%ProcessInterval;
+		if (next_Interval<next)
+			next = next_Interval;
+	}
+
+	//check for FD sample interval
+	if (m_FD_Interval!=0)
+	{
+		int next_Interval = (int)m_FD_Interval - (int)Eng->GetNumberOfTimesteps()%m_FD_Interval;
+		if (next_Interval<next)
+			next = next_Interval;
 	}
-	if (ProcessInterval==0) return next;
-	unsigned int next_Interval = ProcessInterval - Eng->GetNumberOfTimesteps()%ProcessInterval;
-	if (next_Interval<next)
-		next = next_Interval;
 	return next;
 }
 
@@ -87,6 +108,36 @@ void Processing::AddSteps(vector<unsigned int> steps)
 	}
 }
 
+void Processing::AddFrequency(double freq)
+{
+	unsigned int nyquistTS = CalcNyquistNum(freq,Op->GetTimestep());
+
+	if (nyquistTS == 0)
+	{
+		cerr << "Processing::AddFrequency: Requested frequency " << freq << " is too high for the current timestep used... skipping..." << endl;
+		return;
+	}
+	else if (nyquistTS<Op->Exc->GetNyquistNum())
+	{
+		cerr << "Processing::AddFrequency: Warning: Requested frequency " << freq << " is higher than maximum excited frequency..." << endl;
+	}
+
+	if (m_FD_Interval==0)
+		m_FD_Interval = Op->Exc->GetNyquistNum();
+	if (m_FD_Interval>nyquistTS)
+		m_FD_Interval = nyquistTS;
+
+	m_FD_Samples.push_back(freq);
+}
+
+void Processing::AddFrequency(vector<double>  freqs)
+{
+	for (size_t n=0;n<freqs.size();++n)
+	{
+		AddFrequency(freqs.at(n));
+	}
+}
+
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
 	if (Op->SnapToMesh(dstart,start)==false)
@@ -222,11 +273,33 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const
 	file.close();
 }
 
+void Processing::Dump_FD_Data(vector<_Complex double> value, double factor, string filename)
+{
+	ofstream file;
+	file.open( filename.c_str() );
+	if (!file.is_open())
+		cerr << "Can't open file: " << filename << endl;
+
+	for (size_t n=0;n<value.size();++n)
+	{
+		file << m_FD_Samples.at(n) << "\t" << 2.0 * creal(value.at(n))*factor << "\t" << 2.0 * cimag(value.at(n))*factor << "\n";
+	}
+	file.close();
+}
+
 void ProcessingArray::AddProcessing(Processing* proc)
 {
 	ProcessArray.push_back(proc);
 }
 
+void ProcessingArray::FlushNext()
+{
+	for (size_t i=0;i<ProcessArray.size();++i)
+	{
+		ProcessArray.at(i)->FlushNext();
+	}
+}
+
 void ProcessingArray::Reset()
 {
 	for (size_t i=0;i<ProcessArray.size();++i)
diff --git a/FDTD/processing.h b/FDTD/processing.h
index e6a4bf1..ce74068 100644
--- a/FDTD/processing.h
+++ b/FDTD/processing.h
@@ -20,6 +20,8 @@
 
 #include <iostream>
 #include <fstream>
+#include <complex>
+#include <cmath>
 #include "operator.h"
 #include "engine.h"
 
@@ -29,6 +31,7 @@ public:
 	Processing(Operator* op, Engine* eng);
 	virtual ~Processing();
 
+	virtual void Init() {};
 	virtual void Reset();
 
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
@@ -38,6 +41,9 @@ public:
 	void AddStep(unsigned int step);
 	void AddSteps(vector<unsigned int> steps);
 
+	void AddFrequency(double freq);
+	void AddFrequency(vector<double>  freqs);
+
 	bool CheckTimestep();
 	virtual int Process() {return GetNextInterval();}
 
@@ -49,17 +55,24 @@ public:
 	virtual void SetWeight(double weight) {m_weight=weight;}
 	virtual double GetWeight() {return m_weight;}
 
+	//! Invoke this flag to flush all stored data to disk
+	virtual void FlushNext() {m_Flush = true;}
+	virtual void FlushData() {};
+
 	//! Set the dump precision
 	void SetPrecision(unsigned int val) {m_precision = val;}
 
 	virtual void OpenFile( string outfile );
 
 	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
+
 protected:
 	Operator* Op;
 	Engine* Eng;
 	unsigned int m_precision;
 
+	bool m_Flush;
+
 	double m_weight;
 
 	bool Enabled;
@@ -70,6 +83,15 @@ protected:
 	size_t m_PS_pos; //! current position in list of processing steps
 	vector<unsigned int> m_ProcessSteps; //! list of processing steps
 
+	//! Vector of frequency samples
+	vector<double> m_FD_Samples;
+	//! Number of samples already processed
+	unsigned int m_FD_SampleCount;
+	//! Sampling interval needed for the FD_Samples
+	unsigned int m_FD_Interval;
+
+	void Dump_FD_Data(vector<_Complex double> value, double factor, string filename);
+
 	//! define/store snapped start/stop coords as mesh index
 	unsigned int start[3];
 	unsigned int stop[3];
@@ -94,6 +116,9 @@ public:
 
 	void AddProcessing(Processing* proc);
 
+	//! Invoke this flag to flush all stored data to disk for all processings on next Process()
+	void FlushNext();
+
 	void Reset();
 
 	//! Deletes all given processing's, can be helpful, but use carefull!!!
diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp
index ccec3ae..a847f85 100644
--- a/FDTD/processvoltage.cpp
+++ b/FDTD/processvoltage.cpp
@@ -16,6 +16,7 @@
 */
 
 #include "processvoltage.h"
+#include <complex.h>
 #include <iomanip>
 
 ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
@@ -24,16 +25,53 @@ ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
 
 ProcessVoltage::~ProcessVoltage()
 {
+	ProcessVoltage::FlushData();
+}
+
+void ProcessVoltage::Init()
+{
+	FD_voltages.clear();
+	for (size_t n=0;n<m_FD_Samples.size();++n)
+		FD_voltages.push_back(0);
 }
 
 int ProcessVoltage::Process()
 {
 	if (Enabled==false) return -1;
 	if (CheckTimestep()==false) return GetNextInterval();
+
 	FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0);
-//	cerr << voltage << endl;
 	voltage*=m_weight;
-	voltages.push_back(voltage);
-	file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl;
+
+	if (ProcessInterval)
+	{
+		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
+		{
+			voltages.push_back(voltage);
+			file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl;
+		}
+	}
+
+	if (m_FD_Interval)
+	{
+		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0)
+		{
+			double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep();
+			for (size_t n=0;n<m_FD_Samples.size();++n)
+			{
+				FD_voltages.at(n) += voltage * cexp( -2.0 * 1.0i * M_PI * m_FD_Samples.at(n) * T );
+			}
+			++m_FD_SampleCount;
+			if (m_Flush)
+					FlushData();
+			m_Flush = false;
+		}
+	}
+
 	return GetNextInterval();
 }
+
+void ProcessVoltage::FlushData()
+{
+	Dump_FD_Data(FD_voltages,1.0/(double)m_FD_SampleCount,m_filename + "_FD");
+}
diff --git a/FDTD/processvoltage.h b/FDTD/processvoltage.h
index 7272343..7a719dc 100644
--- a/FDTD/processvoltage.h
+++ b/FDTD/processvoltage.h
@@ -27,10 +27,16 @@ public:
 	ProcessVoltage(Operator* op, Engine* eng);
 	virtual ~ProcessVoltage();
 
+	virtual void Init();
+	virtual void FlushData();
+
 	virtual int Process();
 
 protected:
 	vector<FDTD_FLOAT> voltages;
+
+	vector<_Complex double> FD_voltages;
+	void WriteFDVoltages();
 };
 
 #endif // PROCESSVOLTAGE_H
diff --git a/openEMS.pro b/openEMS.pro
index 933c5ea..3df9bad 100644
--- a/openEMS.pro
+++ b/openEMS.pro
@@ -56,7 +56,8 @@ SOURCES += main.cpp \
     FDTD/operator_sse_compressed.cpp \
     FDTD/engine_sse_compressed.cpp \
     FDTD/operator_multithread.cpp \
-    tools/global.cpp
+    tools/global.cpp \
+    tools/useful.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -83,7 +84,8 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/operator_sse_compressed.h \
     FDTD/engine_sse_compressed.h \
     FDTD/operator_multithread.h \
-    tools/global.h
+    tools/global.h \
+    tools/useful.h
 QMAKE_CXXFLAGS_RELEASE = -O3 \
     -g \
 	-march=native
diff --git a/openems.cpp b/openems.cpp
index c210794..1f0293c 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -330,8 +330,10 @@ int openEMS::SetupFDTD(const char* file)
 					continue;
 				}
 				proc->SetProcessInterval(Nyquist/m_OverSampling);
+				proc->AddFrequency(pb->GetFDSamples());
 				proc->DefineStartStopCoord(start,stop);
 				proc->SetWeight(pb->GetWeighting());
+				proc->Init();
 				PA->AddProcessing(proc);
 				prim->SetPrimitiveUsed(true);
 			}
@@ -467,6 +469,8 @@ void openEMS::RunFDTD()
 			cout << " --- Energy: ~" << setw(6) << setprecision(2) << std::scientific << currE << " (decrement: " << setw(6)  << setprecision(2) << std::fixed << fabs(10.0*log10(change)) << "dB)" << endl;
 			prevTime=currTime;
 			prevTS=currTS;
+
+			PA->FlushNext();
 		}
 	}
 
diff --git a/tools/useful.cpp b/tools/useful.cpp
new file mode 100644
index 0000000..8baac9b
--- /dev/null
+++ b/tools/useful.cpp
@@ -0,0 +1,31 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "useful.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <climits>
+
+unsigned int CalcNyquistNum(double fmax, double dT)
+{
+	if (fmax==0) return UINT_MAX;
+	if (dT==0) return 1;
+	double T0 = 1/fmax;
+	return floor(T0/2/dT);
+}
+
diff --git a/tools/useful.h b/tools/useful.h
new file mode 100644
index 0000000..a5a36a3
--- /dev/null
+++ b/tools/useful.h
@@ -0,0 +1,24 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef USEFUL_H
+#define USEFUL_H
+
+//! Calc the nyquist number of timesteps for a given frequency and timestep
+unsigned int CalcNyquistNum(double fmax, double dT);
+
+#endif // USEFUL_H