/*
 * © Copyright 2014 – Urheberrechtshinweis
 * Alle Rechte vorbehalten / All Rights Reserved
 *
 * Programmcode ist urheberrechtlich geschuetzt.
 * Das Urheberrecht liegt, soweit nicht ausdruecklich anders gekennzeichnet, bei Frank Ebner.
 * Keine Verwendung ohne explizite Genehmigung.
 * (vgl. § 106 ff UrhG / § 97 UrhG)
 */

#ifndef FIRREALFACTORY_H
#define FIRREALFACTORY_H

#include "Real.h"

namespace FIR {

	namespace Real {

		class Factory {

			Kernel kernel;

			float sRate_hz;

		public:

			/** ctor */
			Factory(float sRate_hz) : sRate_hz(sRate_hz) {
				;
			}

			/** frequency shift the kernel by multiplying with a frequency */
			void shift(const float shift_hz) {
				for (size_t i = 0; i < kernel.size(); ++i) {
					const float t = (float) i / (float) sRate_hz;
					const float real = std::sin(t * 2 * M_PI * shift_hz);
					kernel[i] = kernel[i] * real;
				}
			}

			/** create a lowpass filte kernel */
			void lowpass(const float cutOff_hz, const int n = 50) {

				kernel.clear();

				for (int i = -n; i <= +n; ++i) {
					const double t = (double) i / (double) sRate_hz;
					const double tmp = 2 * M_PI * cutOff_hz * t;
					const double val = (tmp == 0) ? (1) : (std::sin(tmp) / tmp);
					const double res = val;// * 0.5f;	// why 0.5?
					if (res != res) {throw std::runtime_error("detected NaN");}
					kernel.push_back(res);
				}

			}

			/** apply hamming window to the filter */
			void applyWindowHamming() {
				const int n = (kernel.size()-1)/2;
				int i = -n;
				for (float& f : kernel) {
					f *= winHamming(i+n, n*2);
					++i;
				}
			}

			// https://dsp.stackexchange.com/questions/4693/fir-filter-gain
			/** normalize using the DC-part of the kernel */
			void normalizeDC() {
				float sum = 0;
				for (auto f : kernel) {sum += f;}
				for (auto& f : kernel) {f /= sum;}
			}

			// https://dsp.stackexchange.com/questions/4693/fir-filter-gain
			void normalizeAC(const float freq_hz) {

				const int n = (kernel.size()-1)/2;
				int i = -n;
				float sum = 0;

				for (float f : kernel) {
					const double t = (double) i / (double) sRate_hz;
					const double r = 2 * M_PI * freq_hz * t;
					const double s = std::sin(r);
					sum += f * s;
					++i;
				}

				for (auto& f : kernel) {f /= sum;}

			}





			Kernel getLowpass(const float cutOff_hz, const int n) {
				lowpass(cutOff_hz, n);
				applyWindowHamming();
				normalizeDC();
				return kernel;
			}

			Kernel getBandpass(const float width_hz, const float center_hz, const int n) {
				lowpass(width_hz/2, n);
				applyWindowHamming();
				//normalizeDC();
				shift(center_hz);
				normalizeAC(center_hz);
				return kernel;
			}

			void dumpKernel(const std::string& file, const std::string& varName) {

				std::ofstream out(file);
				out << "# name: " << varName << "\n";
				out << "# type: matrix\n";
				out << "# rows: " << kernel.size() << "\n";
				out << "# columns: 1\n";

				for (const float f : kernel) {
					out << f << "\n";
				}
				out.close();

			}

		private:

			/** get a value from the hamming window */
			static double winHamming(const double t, const double size) {
				return 0.54 - 0.46 * std::cos(2 * M_PI * t / size);
			}

		};

	}

}

#endif // FIRREALFACTORY_H