//////////////////////////////////////////////////////////////////////////
// FTTL.h
//
// Fourier Transform Template Library
// Copyright (C) 2001   Rahul Mittal
// 
// 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 2
// 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, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
//////////////////////////////////////////////////////////////////////////


#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000


#if !defined(__FTTL_H__)
#define __FTTL_H__


#include <vector>
#include <complex>
#include <cmath>
#include <cassert>

#define FTTL_PI		3.14159265358979
#define FTTL_2PI	6.28318530717959


using namespace std;


namespace fttl {

	//////////////////////////////////////////////////////////////////////
	// fourier_analyzer
	//
	template<class T>
		class fourier_analyzer  
	{
		public:
			enum direction {
				reverse = -1,
				forward = 1
			};

			typedef complex<T> complex_t;

		private:
			// Lookup tables
			vector<complex_t>	m_vWpf;	// Forward trigonometric recurrance
			vector<complex_t>	m_vWpr;	// Inverse trigonometric recurrance

			int  bit_reverse(int n, int b) const;
			void compute_trig_recurrance(const int n);
			void transpose(vector<complex_t>& vcData,
						   const int nRows,
						   const int nCols) const;

		public:

			fourier_analyzer(void);
			fourier_analyzer(const direction d);
			fourier_analyzer(const int n);
			virtual ~fourier_analyzer(void);

			//void transform(valarray<complex_t>& acData,
			//			   const direction d);

			void transform(vector<complex_t>& vcData,
						   const direction d);

			void transform(vector<complex_t>& vcData,
						   const int nRows /*rows*/,
						   const int nCols /*cols*/,
						   const direction d);

			void transform(vector<complex_t>::iterator first,
						   vector<complex_t>::iterator last,
						   const direction d);
	};

	template<class T>
		inline int fourier_analyzer<T>::bit_reverse(int n, int b) const
	{
		// Performs a bit reversal of n, starting from the least significant bit
		// and going up to the most significant non-zero bit of b. This
		// bit reversal is key to the Fast Fourier Transform.

		// *** Need to implement a lookup table for bit-reversal to speed things up
		int r = 0;	// return value
		while (b) {
			r <<= 1;
			r |= (n & 1);
			n >>= 1;
			b >>= 1;
		}
		return r;
	}

	template<class T>
		void fourier_analyzer<T>::compute_trig_recurrance(const int n)
	{
		assert ((n & (n - 1)) == 0);	// n must be a power of 2

		int i = n, b = 0;
		while (i) {	// Check if the number of values has changed, to avoid recomputation
			i >>= 1;
			b++;
		}

		if (b > (m_vWpf.size() + 1)) {
		//if (b > m_vWpf.size()) {
			double t, st, ct;
			// Resize to store the # of bits
			m_vWpf.resize(b - 1);
			m_vWpr.resize(b - 1);

			b = 0;
			for (i = 1; i < n; i <<= 1) {
				t = FTTL_PI / i;
				st = sin(t);
				ct = cos(t);
				//m_vWpf[b] = complex_t(cos(t) - 1, sin(t));
				//m_vWpr[b] = complex_t(cos(t) - 1, -sin(t));
				m_vWpf[b] = complex_t(ct - 1, st);
				m_vWpr[b] = complex_t(ct - 1, -st);
				++b;
			}
		}
	}


	template<class T>
		void fourier_analyzer<T>::transpose(vector<complex_t>& vcData,
											const int nRows,
											const int nCols) const
	{
		// Transpose the 1D array that contains 2D data
		assert (vcData.size() == nRows * nCols);

		for (int i = 0; i < nRows; i++) {
			for (int j = 0; j < nCols; j++) {
				swap(vcData[], vcData[]);
			}
		}
	}


	template<class T>
		inline fourier_analyzer<T>::fourier_analyzer(void)
	{
		// Nothing to do
	}

	template<class T>
		inline fourier_analyzer<T>::fourier_analyzer(const int n)
	{
		compute_trig_recurrance(n);
	}


	template<class T>
		inline fourier_analyzer<T>::~fourier_analyzer(void)
	{
		m_vWpf.clear();
		m_vWpr.clear();
	}

	// 1D Complex->Complex Fast Fourier Transform
	template<class T>
		void fourier_analyzer<T>::transform(vector<complex_t>& vcData,
											const direction d)
	{
		const int nSize = vcData.size();
		assert ((nSize & (nSize - 1)) == 0);	// nSize must be a power of 2

		int i, j, k, n, istep, c;
		//double t;
		complex_t wp, w, wt;

		// Step 1: Bit-reversal
		n = nSize >> 1;
		j = 0;
		for (i = 0; i < nSize - 1; i++) {
			j = bit_reverse(i, n);
			if (j > i)
				swap(vcData[j], vcData[i]);
			//k = n;
			//while (k <= j) {
			//	j -= k;
			//	k >>= 1;
			//}
			//j += k;
		}

		compute_trig_recurrance(nSize);

		// Step 2: Danielson-Lanczos formula
		c = 0;
		n = 1;
		while (n < nSize) {
			istep = n << 1;
			//t = d * FTTL_PI / n;
			//wp = complex_t(cos(t) - 1, sin(t));
			wp = (d == direction::forward) ? m_vWpf[c] : m_vWpr[c];
			w = complex_t(1, 0);
			for (k = 0; k < n; k++) {
				for (i = k; i < nSize; i += istep) {
					j = i + n;
					wt = w * vcData[j];
					vcData[j] = vcData[i] - wt;
					vcData[i] += wt;
				}
				w += (w * wp);
			}
			n = istep;
			c++;
		}

		// Scale for the reverse transform
		if (d == direction::reverse) {
			for (vector<complex_t>::iterator i = vcData.begin(); i < vcData.end(); i++)
				*i /= nSize;
		}
	}


	// 1D Complex->Complex Fast Fourier Transform
	template<class T>
		void fourier_analyzer<T>::transform(vector<complex_t>::iterator first,
											vector<complex_t>::iterator last,
											const direction d)
	{
		const int nSize = last - first;
		assert ((nSize & (nSize - 1)) == 0);	// nSize must be a power of 2

		int i, j, k, c, n, istep;
		//double t;
		complex_t wp, w, wt;

		// Step 1: Bit-reversal
		n = nSize >> 1;
		for (i = 0; i < nSize; i++) {
			j = bit_reverse(i, n);
			if (j > i)
				swap(*(first + j), *(first + i));
		}

		compute_trig_recurrance(nSize);

		// Step 2: Danielson-Lanczos formula
		c = 0;
		n = 1;
		while (n < nSize) {
			istep = n << 1;
			//t = d * FTTL_PI / n;
			//wp = complex_t(cos(t) - 1, sin(t));
			wp = (d == direction::forward) ? m_vWpf[c] : m_vWpr[c];
			w = complex_t(1, 0);
			for (k = 0; k < n; k++) {
				for (i = k; i < nSize; i += istep) {
					j = i + n;
					wt = w * *(first + j);
					*(first + j) = *(first + i) - wt;
					*(first + i) += wt;
				}
				w += (w * wp);
			}
			n = istep;
			c++;
		}

		// Scale for the reverse transform
		if (d == direction::reverse) {
			for (vector<complex_t>::iterator i = first; i < last; i++)
				*i /= nSize;
		}
	}


	// 2D Complex->Complex Fast Fourier Transform
	template<class T>
		void fourier_analyzer<T>::transform(vector<complex_t>& vcData,
											const int nRows /*rows*/,
											const int nCols /*cols*/,
											const direction d)
	{
		assert ((nRows & (nRows - 1)) == 0);	// nRows must be a power of 2
		assert ((nCols & (nCols - 1)) == 0);	// nCols must be a power of 2
		assert (vcData.size() == nRows * nCols);

		int i, j, k, row, col, n, ntot, istep, offset, c;
		//double t;
		complex_t wp, w, wt;

		ntot = nRows * nCols;	// == data.size()

		compute_trig_recurrance(nRows > nCols ? nRows : nCols);

		// Divide the algorithm into two stages:
		// 1. Compute fft1d for each row of data
		// 2. Compute fft1d for each column of data (more complicated because of our 2d data is stored as a 1d array)
		// We will not call the fft1d() function because it involves duplicating arrays during Stage 2 which is slow

		// Stage 1 - process rows

		// Step 1: Bit-reversal across each row
		n = nCols >> 1;
		for (row = 0; row < nRows; row++) {
			offset = row * nCols;	// Offset for each row from beginning of array
			j = 0;
			for (col = 0; col < nCols; col++) {
				j = bit_reverse(col, n);
				if (j > col)
					swap(vcData[j + offset], vcData[col + offset]);
				/*
				k = nCols;
				while (k <= j) {
					j -= k;
					k >>= 1;
				}
				j += k;
				*/
			}
		}

		// Step 2: Danielson-Lanczos formula
		for (row = 0; row < nRows; row++) {
			offset = row * nCols;
			c = 0;
			n = 1;
			while (n < nCols) {
				istep = n << 1;
				//t = d * FTTL_PI / n;
				//wp = complex_t(cos(t) - 1, sin(t));
				wp = (d == direction::forward) ? m_vWpf[c] : m_vWpr[c];
				w = complex_t(1, 0);
				for (i = 0; i < n; i++) {
					for (col = i; col < nCols; col += istep) {
						j = col + n;
						wt = w * vcData[j + offset];
						vcData[j + offset] = vcData[col + offset] - wt;
						vcData[col + offset] += wt;
					}
					w += (w * wp);
				}
				n = istep;
				c++;
			}
		}


		// Stage 2 - process columns
		
		// Step 1: Bit-reversal down each column
		n = nRows >> 1;
		for (col = 0; col < nCols; col++) {
			j = 0;
			for (row = 0; row < nRows; row++) {
				j = bit_reverse(row, n);
				if (j > row)
					swap(vcData[j * nCols + col], vcData[row * nCols + col]);
				/*
				k = nRows;
				while (k <= j) {
					j -= k;
					k >>= 1;
				}
				j += k;
				*/
			}
		}

		// Step 2: Danielson-Lanczos formula
		for (col = 0; col < nCols; col++) {
			c = 0;
			n = 1;
			while (n < nRows) {
				istep = n << 1;
				//t = d * FTTL_PI / n;
				//wp = complex_t(cos(t) - 1, sin(t));
				wp = (d == direction::forward) ? m_vWpf[c] : m_vWpr[c];
				w = complex_t(1, 0);
				for (i = 0; i < n; i++) {
					for (row = i; row < nRows; row += istep) {
						j = row + n;
						wt = w * vcData[j * nCols + col];
						vcData[j * nCols + col] = vcData[row * nCols + col] - wt;
						vcData[row * nCols + col] += wt;
					}
					w += (w * wp);
				}
				n = istep;
				c++;
			}
		}

		// Scale for the reverse transform
		if (d == direction::reverse) {
			double d = ntot;
			//d = 1.0;
			for (vector<complex_t>::iterator i = vcData.begin(); i < vcData.end(); i++)
				*i /= d;
		}
	}


	/*
	 * Partial template specialization is unfortunately not supported by MS VC++ 6.0 which
	 * is my only compiler at the moment. I will return to supporting this method when I turn
	 * my linux box on (might take many weeks!) or when MS VC++ 7.0 comes out which I hear
	 * will conform more closely to the standard.
	 */


	/*
	template<class RandomAccessIterator>
		void transform(RandomAccessIterator first, RandomAccessIterator last, const bool isForward)
	{
		iterator_traits<RandomAccessIterator>::distance_type nSize = last - first;
		assert ((nSize & (nSize - 1)) == 0);	// nSize must be a power of 2

		iterator_traits<RandomAccessIterator>::distance_type m, n, iStep;
		iterator_traits<RandomAccessIterator>::value_type wp, w, wt;
		RandomAccessIterator i, j, k;
		double t, pi = isForward ? M_PI : -M_PI;

		// Step 1: Bit-reversal
		n = nSize >> 1;
		for (i = first; i < last; i++) {
			j = bit_reverse<RandomAccessIterator>(i - first, n);
			if (j > i)
				swap(*i, *j);
		}

		// Step 2: Danielson-Lanczos formula
		n = first + 1;
		while (n < nSize) {		// Should execute log_2(nSize) times
			iStep = n * 2;
			t = pi / n;
			wp = iterator_traits<RandomAccessIterator>::value_type(cos(t) - 1, sin(t));
			wp = iterator_traits<RandomAccessIterator>::value_type(1, 0);
			for (k = first; k < n; k++) {
				for (i = k; i < last; i+= iStep) {
					j = i + n;
					wt = w * (*j);
					*j = *i - wt;
					*i += wt;
				}
				w += (w * wp);
			}
			n = iStep;
		}
	}
	*/

	/*
	 * This is a code snippet suggested by [eloi] from #C++ on EFnet which I
	 * am leaving here for my own reference purposes only.
	 */
	/*
	template <typename T> void is_complex (const std::complex<T> &) { }
	
	template <class RandomAccessIterator>
		void fft1d(RandomAccessIterator first, RandomAccessIterator last)
	{
		is_complex(std::iterator_traits<RandomAccessIterator>::value_type());
	}
	*/

	// Code snippets from #C++ on EFNet
	// template <class T> std::string ToString ( T &t ) { ostringstream os; os << t; return os.str(); }


}	// namespace fttl


#endif	// !defined(__FTTL_H__)