BogaudioModules

BogaudioModules for VCV Rack
Log | Files | Refs | README | LICENSE

FFTRealPassDirect.hpp (5654B)


      1 /*****************************************************************************
      2 
      3         FFTRealPassDirect.hpp
      4         By Laurent de Soras
      5 
      6 --- Legal stuff ---
      7 
      8 This program is free software. It comes without any warranty, to
      9 the extent permitted by applicable law. You can redistribute it
     10 and/or modify it under the terms of the Do What The Fuck You Want
     11 To Public License, Version 2, as published by Sam Hocevar. See
     12 http://sam.zoy.org/wtfpl/COPYING for more details.
     13 
     14 *Tab=3***********************************************************************/
     15 
     16 
     17 
     18 #if defined (ffft_FFTRealPassDirect_CURRENT_CODEHEADER)
     19 	#error Recursive inclusion of FFTRealPassDirect code header.
     20 #endif
     21 #define	ffft_FFTRealPassDirect_CURRENT_CODEHEADER
     22 
     23 #if ! defined (ffft_FFTRealPassDirect_CODEHEADER_INCLUDED)
     24 #define	ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
     25 
     26 
     27 
     28 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
     29 
     30 #include	"ffft/FFTRealUseTrigo.h"
     31 
     32 
     33 
     34 namespace ffft
     35 {
     36 
     37 
     38 
     39 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
     40 
     41 
     42 
     43 template <>
     44 inline void	FFTRealPassDirect <1>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
     45 {
     46 	// First and second pass at once
     47 	const long		qlen = len >> 2;
     48 
     49 	long				coef_index = 0;
     50 	do
     51 	{
     52 		// To do: unroll the loop (2x).
     53 		const long		ri_0 = br_ptr [coef_index >> 2];
     54 		const long		ri_1 = ri_0 + 2 * qlen;	// bit_rev_lut_ptr [coef_index + 1];
     55 		const long		ri_2 = ri_0 + 1 * qlen;	// bit_rev_lut_ptr [coef_index + 2];
     56 		const long		ri_3 = ri_0 + 3 * qlen;	// bit_rev_lut_ptr [coef_index + 3];
     57 
     58 		DataType	* const	df2 = dest_ptr + coef_index;
     59 		df2 [1] = x_ptr [ri_0] - x_ptr [ri_1];
     60 		df2 [3] = x_ptr [ri_2] - x_ptr [ri_3];
     61 
     62 		const DataType	sf_0 = x_ptr [ri_0] + x_ptr [ri_1];
     63 		const DataType	sf_2 = x_ptr [ri_2] + x_ptr [ri_3];
     64 
     65 		df2 [0] = sf_0 + sf_2;
     66 		df2 [2] = sf_0 - sf_2;
     67 
     68 		coef_index += 4;
     69 	}
     70 	while (coef_index < len);
     71 }
     72 
     73 template <>
     74 inline void	FFTRealPassDirect <2>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
     75 {
     76 	// Executes "previous" passes first. Inverts source and destination buffers
     77 	FFTRealPassDirect <1>::process (
     78 		len,
     79 		src_ptr,
     80 		dest_ptr,
     81 		x_ptr,
     82 		cos_ptr,
     83 		cos_len,
     84 		br_ptr,
     85 		osc_list
     86 	);
     87 
     88 	// Third pass
     89 	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
     90 
     91 	long				coef_index = 0;
     92 	do
     93 	{
     94 		dest_ptr [coef_index    ] = src_ptr [coef_index] + src_ptr [coef_index + 4];
     95 		dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];
     96 		dest_ptr [coef_index + 2] = src_ptr [coef_index + 2];
     97 		dest_ptr [coef_index + 6] = src_ptr [coef_index + 6];
     98 
     99 		DataType			v;
    100 
    101 		v = (src_ptr [coef_index + 5] - src_ptr [coef_index + 7]) * sqrt2_2;
    102 		dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + v;
    103 		dest_ptr [coef_index + 3] = src_ptr [coef_index + 1] - v;
    104 
    105 		v = (src_ptr [coef_index + 5] + src_ptr [coef_index + 7]) * sqrt2_2;
    106 		dest_ptr [coef_index + 5] = v + src_ptr [coef_index + 3];
    107 		dest_ptr [coef_index + 7] = v - src_ptr [coef_index + 3];
    108 
    109 		coef_index += 8;
    110 	}
    111 	while (coef_index < len);
    112 }
    113 
    114 template <int PASS>
    115 void	FFTRealPassDirect <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])
    116 {
    117 	// Executes "previous" passes first. Inverts source and destination buffers
    118 	FFTRealPassDirect <PASS - 1>::process (
    119 		len,
    120 		src_ptr,
    121 		dest_ptr,
    122 		x_ptr,
    123 		cos_ptr,
    124 		cos_len,
    125 		br_ptr,
    126 		osc_list
    127 	);
    128 
    129 	const long		dist = 1L << (PASS - 1);
    130 	const long		c1_r = 0;
    131 	const long		c1_i = dist;
    132 	const long		c2_r = dist * 2;
    133 	const long		c2_i = dist * 3;
    134 	const long		cend = dist * 4;
    135 	const long		table_step = cos_len >> (PASS - 1);
    136 
    137    enum {	TRIGO_OSC		= PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT	};
    138 	enum {	TRIGO_DIRECT	= (TRIGO_OSC >= 0) ? 1 : 0	};
    139 
    140 	long				coef_index = 0;
    141 	do
    142 	{
    143 		const DataType	* const	sf = src_ptr + coef_index;
    144 		DataType			* const	df = dest_ptr + coef_index;
    145 
    146 		// Extreme coefficients are always real
    147 		df [c1_r] = sf [c1_r] + sf [c2_r];
    148 		df [c2_r] = sf [c1_r] - sf [c2_r];
    149 		df [c1_i] = sf [c1_i];
    150 		df [c2_i] = sf [c2_i];
    151 
    152 		FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);
    153 
    154 		// Others are conjugate complex numbers
    155 		for (long i = 1; i < dist; ++ i)
    156 		{
    157 			DataType			c;
    158 			DataType			s;
    159 			FFTRealUseTrigo <TRIGO_DIRECT>::iterate (
    160 				osc_list [TRIGO_OSC],
    161 				c,
    162 				s,
    163 				cos_ptr,
    164 				i * table_step,
    165 				(dist - i) * table_step
    166 			);
    167 
    168 			const DataType	sf_r_i = sf [c1_r + i];
    169 			const DataType	sf_i_i = sf [c1_i + i];
    170 
    171 			const DataType	v1 = sf [c2_r + i] * c - sf [c2_i + i] * s;
    172 			df [c1_r + i] = sf_r_i + v1;
    173 			df [c2_r - i] = sf_r_i - v1;
    174 
    175 			const DataType	v2 = sf [c2_r + i] * s + sf [c2_i + i] * c;
    176 			df [c2_r + i] = v2 + sf_i_i;
    177 			df [cend - i] = v2 - sf_i_i;
    178 		}
    179 
    180 		coef_index += cend;
    181 	}
    182 	while (coef_index < len);
    183 }
    184 
    185 
    186 
    187 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
    188 
    189 
    190 
    191 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
    192 
    193 
    194 
    195 }	// namespace ffft
    196 
    197 
    198 
    199 #endif	// ffft_FFTRealPassDirect_CODEHEADER_INCLUDED
    200 
    201 #undef ffft_FFTRealPassDirect_CURRENT_CODEHEADER
    202 
    203 
    204 
    205 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/