/*****************************************************************************

        TestAccuracy.hpp

        Copyright (c) 2005 Laurent de Soras

--- Legal stuff ---

This library is free software; you can redistribute it and/or

modify it under the terms of the GNU Lesser General Public

License as published by the Free Software Foundation; either

version 2.1 of the License, or (at your option) any later version.

This library 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

Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public

License along with this library; if not, write to the Free Software

Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

*Tab=3***********************************************************************/

#if defined (TestAccuracy_CURRENT_CODEHEADER)

	#error Recursive inclusion of TestAccuracy code header.

#endif

#define	TestAccuracy_CURRENT_CODEHEADER

#if ! defined (TestAccuracy_CODEHEADER_INCLUDED)

#define	TestAccuracy_CODEHEADER_INCLUDED

/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

#include	"def.h"

#include	"test_fnc.h"

#include	"TestWhiteNoiseGen.h"

#include	<typeinfo>

#include	<vector>

#include	<cmath>

#include	<cstdio>

static const double	TestAccuracy_LN10	= 2.3025850929940456840179914546844;

/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

template <class FO>

int	TestAccuracy <FO>::perform_test_single_object (FO &fft)

{

	assert (&fft != 0);

	using namespace std;

	int				ret_val = 0;

	const std::type_info &	ti = typeid (fft);

	const char *	class_name_0 = ti.name ();

	if (ret_val == 0)

	{

		ret_val = perform_test_d (fft, class_name_0);

	}

	if (ret_val == 0)

	{

		ret_val = perform_test_i (fft, class_name_0);

	}

	if (ret_val == 0)

	{

		ret_val = perform_test_di (fft, class_name_0);

	}

	if (ret_val == 0)

	{

		printf ("\n");

	}

	return (ret_val);

}

template <class FO>

int	TestAccuracy <FO>::perform_test_d (FO &fft, const char *class_name_0)

{

	assert (&fft != 0);

   assert (class_name_0 != 0);

	using namespace std;

	int				ret_val = 0;

	const long		len = fft.get_length ();

   const long     nbr_tests = limit (

      NBR_ACC_TESTS / len / len,

      1L,

      static_cast <long> (MAX_NBR_TESTS)

   );

	printf ("Testing %s::do_fft () [%ld samples]... ", class_name_0, len);

	fflush (stdout);

	TestWhiteNoiseGen <DataType>	noise;

	std::vector <DataType>	x (len);

	std::vector <DataType>	s1 (len);

	std::vector <DataType>	s2 (len);

	BigFloat			err_avg = 0;

	for (long test = 0; test < nbr_tests && ret_val == 0; ++ test)

	{

		noise.generate (&x [0], len);

		fft.do_fft (&s1 [0], &x [0]);

		compute_tf (&s2 [0], &x [0], len);

		BigFloat			max_err;

		compare_vect_display (&s1 [0], &s2 [0], len, max_err);

		err_avg += max_err;

	}

	err_avg /= NBR_ACC_TESTS;

	printf ("done.\n");

	printf (

		"Average maximum error: %.6f %% (%f dB)\n",

		static_cast <double> (err_avg * 100),

		static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg))

	);

	return (ret_val);

}

template <class FO>

int	TestAccuracy <FO>::perform_test_i (FO &fft, const char *class_name_0)

{

	assert (&fft != 0);

   assert (class_name_0 != 0);

	using namespace std;

	int				ret_val = 0;

	const long		len = fft.get_length ();

   const long     nbr_tests = limit (

      NBR_ACC_TESTS / len / len,

      10L,

      static_cast <long> (MAX_NBR_TESTS)

   );

	printf ("Testing %s::do_ifft () [%ld samples]... ", class_name_0, len);

	fflush (stdout);

	TestWhiteNoiseGen <DataType>	noise;

	std::vector <DataType>	s (len);

	std::vector <DataType>	x1 (len);

	std::vector <DataType>	x2 (len);

	BigFloat			err_avg = 0;

	for (long test = 0; test < nbr_tests && ret_val == 0; ++ test)

	{

		noise.generate (&s [0], len);

		fft.do_ifft (&s [0], &x1 [0]);

		compute_itf (&x2 [0], &s [0], len);

		BigFloat			max_err;

		compare_vect_display (&x1 [0], &x2 [0], len, max_err);

		err_avg += max_err;

	}

	err_avg /= NBR_ACC_TESTS;

	printf ("done.\n");

	printf (

		"Average maximum error: %.6f %% (%f dB)\n",

		static_cast <double> (err_avg * 100),

		static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg))

	);

	return (ret_val);

}

template <class FO>

int	TestAccuracy <FO>::perform_test_di (FO &fft, const char *class_name_0)

{

	assert (&fft != 0);

   assert (class_name_0 != 0);

	using namespace std;

	int				ret_val = 0;

	const long		len = fft.get_length ();

   const long     nbr_tests = limit (

      NBR_ACC_TESTS / len / len,

      1L,

      static_cast <long> (MAX_NBR_TESTS)

   );

	printf (

		"Testing %s::do_fft () / do_ifft () / rescale () [%ld samples]... ",

		class_name_0,

		len

	);

	fflush (stdout);

	TestWhiteNoiseGen <DataType>	noise;

	std::vector <DataType>	x (len);

	std::vector <DataType>	s (len);

	std::vector <DataType>	y (len);

	BigFloat			err_avg = 0;

	for (long test = 0; test < nbr_tests && ret_val == 0; ++ test)

	{

		noise.generate (&x [0], len);

		fft.do_fft (&s [0], &x [0]);

		fft.do_ifft (&s [0], &y [0]);

		fft.rescale (&y [0]);

		BigFloat			max_err;

		compare_vect_display (&x [0], &y [0], len, max_err);

		err_avg += max_err;

	}

	err_avg /= NBR_ACC_TESTS;

	printf ("done.\n");

	printf (

		"Average maximum error: %.6f %% (%f dB)\n",

		static_cast <double> (err_avg * 100),

		static_cast <double> ((20 / TestAccuracy_LN10) * log (err_avg))

	);

	return (ret_val);

}

/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

// Positive transform

template <class FO>

void	TestAccuracy <FO>::compute_tf (DataType s [], const DataType x [], long length)

{

	assert (s != 0);

	assert (x != 0);

	assert (length >= 2);

	assert ((length & 1) == 0);

	const long		nbr_bins = length >> 1;

	// DC and Nyquist

	BigFloat			dc = 0;

	BigFloat			ny = 0;

	for (long pos = 0; pos < length; pos += 2)

	{

		const BigFloat	even = x [pos    ];

		const BigFloat	odd  = x [pos + 1];

		dc += even + odd;

		ny += even - odd;

	}

	s [0       ] = static_cast <DataType> (dc);

	s [nbr_bins] = static_cast <DataType> (ny);

	// Regular bins

	for (long bin = 1; bin < nbr_bins; ++ bin)

	{

		BigFloat			sum_r = 0;

		BigFloat			sum_i = 0;

		const BigFloat	m = bin * static_cast <BigFloat> (2 * PI) / length;

		for (long pos = 0; pos < length; ++pos)

		{

			using namespace std;

			const BigFloat	phase = pos * m;

			const BigFloat	e_r = cos (phase);

			const BigFloat	e_i = sin (phase);

			sum_r += x [pos] * e_r;

			sum_i += x [pos] * e_i;

		}

		s [           bin] = static_cast <DataType> (sum_r);

		s [nbr_bins + bin] = static_cast <DataType> (sum_i);

	}

}

// Negative transform

template <class FO>

void	TestAccuracy <FO>::compute_itf (DataType x [], const DataType s [], long length)

{

	assert (s != 0);

	assert (x != 0);

	assert (length >= 2);

	assert ((length & 1) == 0);

	const long		nbr_bins = length >> 1;

	// DC and Nyquist

	BigFloat			dc = s [0       ];

	BigFloat			ny = s [nbr_bins];

	// Regular bins

	for (long pos = 0; pos < length; ++pos)

	{

		BigFloat				sum = dc + ny * (1 - 2 * (pos & 1));

		const BigFloat		m = pos * static_cast <BigFloat> (-2 * PI) / length;

		for (long bin = 1; bin < nbr_bins; ++ bin)

		{

			using namespace std;

			const BigFloat	phase = bin * m;

			const BigFloat	e_r = cos (phase);

			const BigFloat	e_i = sin (phase);

			sum += 2 * (  e_r * s [bin           ]

			            - e_i * s [bin + nbr_bins]);

		}

		x [pos] = static_cast <DataType> (sum);

	}

}

template <class FO>

int	TestAccuracy <FO>::compare_vect_display (const DataType x_ptr [], const DataType y_ptr [], long len, BigFloat &max_err_rel)

{

	assert (x_ptr != 0);

	assert (y_ptr != 0);

	assert (len > 0);

	assert (&max_err_rel != 0);

	using namespace std;

	int				ret_val = 0;

	BigFloat			power = compute_power (&x_ptr [0], &y_ptr [0], len);

	BigFloat			power_dif;

	long				max_err_pos;

	compare_vect (&x_ptr [0], &y_ptr [0], power_dif, max_err_pos, len);

	if (power == 0)

	{

		power = power_dif;

	}

	const BigFloat	power_err_rel = power_dif / power;

	BigFloat	      max_err = 0;

	max_err_rel = 0;

	if (max_err_pos >= 0)

	{

      max_err = y_ptr [max_err_pos] - x_ptr [max_err_pos];

		max_err_rel = 2 * fabs (max_err) / (  fabs (y_ptr [max_err_pos])

		                                    + fabs (x_ptr [max_err_pos]));

	}

	if (power_err_rel > 0.001)

	{

		printf ("Power error  : %f (%.6f %%)\n",

			static_cast <double> (power_err_rel),

			static_cast <double> (power_err_rel * 100)

		);

		if (max_err_pos >= 0)

		{

			printf (

				"Maximum error: %f - %f = %f (%f)\n",

				static_cast <double> (y_ptr [max_err_pos]),

				static_cast <double> (x_ptr [max_err_pos]),

				static_cast <double> (max_err),

				static_cast <double> (max_err_pos)

			);

		}

	}

	return (ret_val);

}

template <class FO>

typename TestAccuracy <FO>::BigFloat	TestAccuracy <FO>::compute_power (const DataType x_ptr [], long len)

{

	assert (x_ptr != 0);

	assert (len > 0);

	BigFloat		power = 0;

	for (long pos = 0; pos < len; ++pos)

	{

		const BigFloat	val = x_ptr [pos];

		power += val * val;

	}

	using namespace std;

	power = sqrt (power) / len;

	return (power);

}

template <class FO>

typename TestAccuracy <FO>::BigFloat	TestAccuracy <FO>::compute_power (const DataType x_ptr [], const DataType y_ptr [], long len)

{

	assert (x_ptr != 0);

	assert (y_ptr != 0);

	assert (len > 0);

	return ((compute_power (x_ptr, len) + compute_power (y_ptr, len)) * 0.5);

}

template <class FO>

void	TestAccuracy <FO>::compare_vect (const DataType x_ptr [], const DataType y_ptr [], BigFloat &power, long &max_err_pos, long len)

{

	assert (x_ptr != 0);

	assert (y_ptr != 0);

	assert (len > 0);

	assert (&power != 0);

	assert (&max_err_pos != 0);

	power = 0;

	BigFloat			max_dif2 = 0;

	max_err_pos = -1;

	for (long pos = 0; pos < len; ++pos)

	{

		const BigFloat	x = x_ptr [pos];

		const BigFloat	y = y_ptr [pos];

		const BigFloat	dif = y - x;

		const BigFloat	dif2 = dif * dif;

		power += dif2;

		if (dif2 > max_dif2)

		{

			max_err_pos = pos;

			max_dif2 = dif2;

		}

	}

	using namespace std;

	power = sqrt (power) / len;

}

#endif	// TestAccuracy_CODEHEADER_INCLUDED

#undef TestAccuracy_CURRENT_CODEHEADER

/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/