demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp
changeset 25 e24348a560a6
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp	Fri Jun 11 14:24:45 2010 +0300
@@ -0,0 +1,229 @@
+/*****************************************************************************
+
+        FFTRealPassInverse.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 (FFTRealPassInverse_CURRENT_CODEHEADER)
+	#error Recursive inclusion of FFTRealPassInverse code header.
+#endif
+#define	FFTRealPassInverse_CURRENT_CODEHEADER
+
+#if ! defined (FFTRealPassInverse_CODEHEADER_INCLUDED)
+#define	FFTRealPassInverse_CODEHEADER_INCLUDED
+
+
+
+/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
+
+#include	"FFTRealUseTrigo.h"
+
+
+
+/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
+
+
+
+template <int PASS>
+void	FFTRealPassInverse <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType f_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	process_internal (
+		len,
+		dest_ptr,
+		f_ptr,
+		cos_ptr,
+		cos_len,
+		br_ptr,
+		osc_list
+	);
+	FFTRealPassInverse <PASS - 1>::process_rec (
+		len,
+		src_ptr,
+		dest_ptr,
+		cos_ptr,
+		cos_len,
+		br_ptr,
+		osc_list
+	);
+}
+
+
+
+template <int PASS>
+void	FFTRealPassInverse <PASS>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	process_internal (
+		len,
+		dest_ptr,
+		src_ptr,
+		cos_ptr,
+		cos_len,
+		br_ptr,
+		osc_list
+	);
+	FFTRealPassInverse <PASS - 1>::process_rec (
+		len,
+		src_ptr,
+		dest_ptr,
+		cos_ptr,
+		cos_len,
+		br_ptr,
+		osc_list
+	);
+}
+
+template <>
+void	FFTRealPassInverse <0>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	// Stops recursion
+}
+
+
+
+template <int PASS>
+void	FFTRealPassInverse <PASS>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	const long		dist = 1L << (PASS - 1);
+	const long		c1_r = 0;
+	const long		c1_i = dist;
+	const long		c2_r = dist * 2;
+	const long		c2_i = dist * 3;
+	const long		cend = dist * 4;
+	const long		table_step = cos_len >> (PASS - 1);
+
+   enum {	TRIGO_OSC		= PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT	};
+	enum {	TRIGO_DIRECT	= (TRIGO_OSC >= 0) ? 1 : 0	};
+
+	long				coef_index = 0;
+	do
+	{
+		const DataType	* const	sf = src_ptr + coef_index;
+		DataType			* const	df = dest_ptr + coef_index;
+
+		// Extreme coefficients are always real
+		df [c1_r] = sf [c1_r] + sf [c2_r];
+		df [c2_r] = sf [c1_r] - sf [c2_r];
+		df [c1_i] = sf [c1_i] * 2;
+		df [c2_i] = sf [c2_i] * 2;
+
+		FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
+
+		// Others are conjugate complex numbers
+		for (long i = 1; i < dist; ++ i)
+		{
+			df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i];
+			df [c1_i + i] = sf [c2_r + i] - sf [cend - i];
+
+			DataType			c;
+			DataType			s;
+			FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
+				osc_list [TRIGO_OSC],
+				c,
+				s,
+				cos_ptr,
+				i * table_step,
+				(dist - i) * table_step
+			);
+
+			const DataType	vr = sf [c1_r + i] - sf [c2_r - i];
+			const DataType	vi = sf [c2_r + i] + sf [cend - i];
+
+			df [c2_r + i] = vr * c + vi * s;
+			df [c2_i + i] = vi * c - vr * s;
+		}
+
+		coef_index += cend;
+	}
+	while (coef_index < len);
+}
+
+template <>
+void	FFTRealPassInverse <2>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	// Antepenultimate pass
+	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
+
+	long				coef_index = 0;
+	do
+	{
+		dest_ptr [coef_index    ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
+		dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
+		dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2;
+		dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2;
+
+		dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3];
+		dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7];
+
+		const DataType	vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3];
+		const DataType	vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7];
+
+		dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2;
+		dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2;
+
+		coef_index += 8;
+	}
+	while (coef_index < len);
+}
+
+template <>
+void	FFTRealPassInverse <1>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
+{
+	// Penultimate and last pass at once
+	const long		qlen = len >> 2;
+
+	long				coef_index = 0;
+	do
+	{
+		const long		ri_0 = br_ptr [coef_index >> 2];
+
+		const DataType	b_0 = src_ptr [coef_index    ] + src_ptr [coef_index + 2];
+		const DataType	b_2 = src_ptr [coef_index    ] - src_ptr [coef_index + 2];
+		const DataType	b_1 = src_ptr [coef_index + 1] * 2;
+		const DataType	b_3 = src_ptr [coef_index + 3] * 2;
+
+		dest_ptr [ri_0           ] = b_0 + b_1;
+		dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1;
+		dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3;
+		dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3;
+
+		coef_index += 4;
+	}
+	while (coef_index < len);
+}
+
+
+
+/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
+
+
+
+/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
+
+
+
+#endif	// FFTRealPassInverse_CODEHEADER_INCLUDED
+
+#undef FFTRealPassInverse_CURRENT_CODEHEADER
+
+
+
+/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/