CDSPHBDownsampler.h 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401
  1. //$ nobt
  2. //$ nocpp
  3. /**
  4. * @file CDSPHBDownsampler.h
  5. *
  6. * @brief Half-band downsampling convolver class.
  7. *
  8. * This file includes half-band downsampling convolver class.
  9. *
  10. * r8brain-free-src Copyright (c) 2013-2022 Aleksey Vaneev
  11. * See the "LICENSE" file for license.
  12. */
  13. #ifndef R8B_CDSPHBDOWNSAMPLER_INCLUDED
  14. #define R8B_CDSPHBDOWNSAMPLER_INCLUDED
  15. #include "CDSPHBUpsampler.h"
  16. namespace r8b {
  17. /**
  18. * @brief Half-band downsampler class.
  19. *
  20. * Class implements brute-force half-band 2X downsampling that uses small
  21. * sparse symmetric FIR filters. The output has 2.0 gain.
  22. */
  23. class CDSPHBDownsampler : public CDSPProcessor
  24. {
  25. public:
  26. /**
  27. * Constructor initalizes the half-band downsampler.
  28. *
  29. * @param ReqAtten Required half-band filter attentuation.
  30. * @param SteepIndex Steepness index - 0=steepest. Corresponds to general
  31. * downsampling ratio, e.g. at 4x downsampling 0 is used, at 8x
  32. * downsampling 1 is used, etc.
  33. * @param IsThird "True" if 1/3 resampling is performed.
  34. * @param PrevLatency Latency, in samples (any value >=0), which was left
  35. * in the output signal by a previous process. Whole-number latency will
  36. * be consumed by *this object while remaining fractional latency can be
  37. * obtained via the getLatencyFrac() function.
  38. */
  39. CDSPHBDownsampler( const double ReqAtten, const int SteepIndex,
  40. const bool IsThird, const double PrevLatency )
  41. {
  42. static const CConvolveFn FltConvFn[ 14 ] = {
  43. &CDSPHBDownsampler :: convolve1, &CDSPHBDownsampler :: convolve2,
  44. &CDSPHBDownsampler :: convolve3, &CDSPHBDownsampler :: convolve4,
  45. &CDSPHBDownsampler :: convolve5, &CDSPHBDownsampler :: convolve6,
  46. &CDSPHBDownsampler :: convolve7, &CDSPHBDownsampler :: convolve8,
  47. &CDSPHBDownsampler :: convolve9, &CDSPHBDownsampler :: convolve10,
  48. &CDSPHBDownsampler :: convolve11, &CDSPHBDownsampler :: convolve12,
  49. &CDSPHBDownsampler :: convolve13,
  50. &CDSPHBDownsampler :: convolve14 };
  51. int fltt;
  52. double att;
  53. if( IsThird )
  54. {
  55. CDSPHBUpsampler :: getHBFilterThird( ReqAtten, SteepIndex, fltp,
  56. fltt, att );
  57. }
  58. else
  59. {
  60. CDSPHBUpsampler :: getHBFilter( ReqAtten, SteepIndex, fltp, fltt,
  61. att );
  62. }
  63. convfn = FltConvFn[ fltt - 1 ];
  64. fll = fltt * 2 - 1;
  65. fl2 = fll;
  66. flo = fll + fl2;
  67. flb = BufLen - fll;
  68. BufRP = Buf + fll;
  69. LatencyFrac = PrevLatency * 0.5;
  70. Latency = (int) LatencyFrac;
  71. LatencyFrac -= Latency;
  72. R8BASSERT( Latency >= 0 );
  73. R8BCONSOLE( "CDSPHBDownsampler: taps=%i third=%i att=%.1f io=1/2\n",
  74. fltt, (int) IsThird, att );
  75. clear();
  76. }
  77. virtual int getLatency() const
  78. {
  79. return( 0 );
  80. }
  81. virtual double getLatencyFrac() const
  82. {
  83. return( LatencyFrac );
  84. }
  85. virtual int getMaxOutLen( const int MaxInLen ) const
  86. {
  87. R8BASSERT( MaxInLen >= 0 );
  88. return(( MaxInLen + 1 ) >> 1 );
  89. }
  90. virtual void clear()
  91. {
  92. LatencyLeft = Latency;
  93. BufLeft = 0;
  94. WritePos = 0;
  95. ReadPos = flb; // Set "read" position to account for filter's latency.
  96. memset( &Buf[ ReadPos ], 0, ( BufLen - flb ) * sizeof( Buf[ 0 ]));
  97. }
  98. virtual int process( double* ip, int l, double*& op0 )
  99. {
  100. R8BASSERT( l >= 0 );
  101. double* op = op0;
  102. while( l > 0 )
  103. {
  104. // Add new input samples to both halves of the ring buffer.
  105. const int b = min( min( l, BufLen - WritePos ), flb - BufLeft );
  106. double* const wp1 = Buf + WritePos;
  107. memcpy( wp1, ip, b * sizeof( wp1[ 0 ]));
  108. if( WritePos < flo )
  109. {
  110. const int c = min( b, flo - WritePos );
  111. memcpy( wp1 + BufLen, wp1, c * sizeof( wp1[ 0 ]));
  112. }
  113. ip += b;
  114. WritePos = ( WritePos + b ) & BufLenMask;
  115. l -= b;
  116. BufLeft += b;
  117. // Produce output.
  118. if( BufLeft > fl2 )
  119. {
  120. const int c = ( BufLeft - fl2 + 1 ) >> 1;
  121. double* const opend = op + c;
  122. ( *convfn )( op, opend, fltp, BufRP, ReadPos );
  123. op = opend;
  124. const int c2 = c + c;
  125. ReadPos = ( ReadPos + c2 ) & BufLenMask;
  126. BufLeft -= c2;
  127. }
  128. }
  129. int ol = (int) ( op - op0 );
  130. if( LatencyLeft != 0 )
  131. {
  132. if( LatencyLeft >= ol )
  133. {
  134. LatencyLeft -= ol;
  135. return( 0 );
  136. }
  137. ol -= LatencyLeft;
  138. op0 += LatencyLeft;
  139. LatencyLeft = 0;
  140. }
  141. return( ol );
  142. }
  143. private:
  144. static const int BufLenBits = 10; ///< The length of the ring buffer,
  145. ///< expressed as Nth power of 2. This value can be reduced if it is
  146. ///< known that only short input buffers will be passed to the
  147. ///< interpolator. The minimum value of this parameter is 5, and
  148. ///< 1 << BufLenBits should be at least 3 times larger than the
  149. ///< FilterLen.
  150. ///<
  151. static const int BufLen = 1 << BufLenBits; ///< The length of the ring
  152. ///< buffer. The actual length is twice as long to allow "beyond max
  153. ///< position" positioning.
  154. ///<
  155. static const int BufLenMask = BufLen - 1; ///< Mask used for quick buffer
  156. ///< position wrapping.
  157. ///<
  158. double Buf[ BufLen + 54 ]; ///< The ring buffer, including overrun
  159. ///< protection for the largest filter.
  160. ///<
  161. const double* fltp; ///< Half-band filter taps.
  162. ///<
  163. int fll; ///< Input latency.
  164. ///<
  165. int fl2; ///< Right-side filter length.
  166. ///<
  167. int flo; ///< Overrun length.
  168. ///<
  169. int flb; ///< Initial read position and maximal buffer write length.
  170. ///<
  171. const double* BufRP; ///< Offseted Buf pointer at ReadPos=0.
  172. ///<
  173. int Latency; ///< Initial latency that should be removed from the output.
  174. ///<
  175. double LatencyFrac; ///< Fractional latency left on the output.
  176. ///<
  177. int BufLeft; ///< The number of samples left in the buffer to process.
  178. ///< When this value is below FilterLenD2Plus1, the interpolation
  179. ///< cycle ends.
  180. ///<
  181. int WritePos; ///< The current buffer write position. Incremented together
  182. ///< with the BufLeft variable.
  183. ///<
  184. int ReadPos; ///< The current buffer read position.
  185. ///<
  186. int LatencyLeft; ///< Latency left to remove.
  187. ///<
  188. typedef void( *CConvolveFn )( double* op, double* const opend,
  189. const double* const flt, const double* const rp0, int rpos ); ///<
  190. ///< Convolution function type.
  191. ///<
  192. CConvolveFn convfn; ///< Convolution function in use.
  193. ///<
  194. #define R8BHBC1( fn ) \
  195. static void fn( double* op, double* const opend, const double* const flt, \
  196. const double* const rp0, int rpos ) \
  197. { \
  198. while( op != opend ) \
  199. { \
  200. const double* const rp = rp0 + rpos; \
  201. *op = rp[ 0 ] +
  202. #define R8BHBC2 \
  203. rpos = ( rpos + 2 ) & BufLenMask; \
  204. op++; \
  205. } \
  206. }
  207. R8BHBC1( convolve1 )
  208. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]);
  209. R8BHBC2
  210. R8BHBC1( convolve2 )
  211. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  212. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]);
  213. R8BHBC2
  214. R8BHBC1( convolve3 )
  215. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  216. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  217. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]);
  218. R8BHBC2
  219. R8BHBC1( convolve4 )
  220. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  221. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  222. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  223. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]);
  224. R8BHBC2
  225. R8BHBC1( convolve5 )
  226. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  227. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  228. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  229. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  230. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]);
  231. R8BHBC2
  232. R8BHBC1( convolve6 )
  233. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  234. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  235. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  236. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  237. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  238. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]);
  239. R8BHBC2
  240. R8BHBC1( convolve7 )
  241. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  242. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  243. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  244. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  245. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  246. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  247. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]);
  248. R8BHBC2
  249. R8BHBC1( convolve8 )
  250. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  251. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  252. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  253. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  254. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  255. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  256. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  257. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]);
  258. R8BHBC2
  259. R8BHBC1( convolve9 )
  260. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  261. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  262. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  263. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  264. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  265. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  266. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  267. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  268. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]);
  269. R8BHBC2
  270. R8BHBC1( convolve10 )
  271. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  272. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  273. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  274. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  275. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  276. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  277. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  278. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  279. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]) +
  280. flt[ 9 ] * ( rp[ 19 ] + rp[ -19 ]);
  281. R8BHBC2
  282. R8BHBC1( convolve11 )
  283. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  284. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  285. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  286. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  287. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  288. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  289. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  290. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  291. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]) +
  292. flt[ 9 ] * ( rp[ 19 ] + rp[ -19 ]) +
  293. flt[ 10 ] * ( rp[ 21 ] + rp[ -21 ]);
  294. R8BHBC2
  295. R8BHBC1( convolve12 )
  296. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  297. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  298. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  299. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  300. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  301. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  302. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  303. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  304. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]) +
  305. flt[ 9 ] * ( rp[ 19 ] + rp[ -19 ]) +
  306. flt[ 10 ] * ( rp[ 21 ] + rp[ -21 ]) +
  307. flt[ 11 ] * ( rp[ 23 ] + rp[ -23 ]);
  308. R8BHBC2
  309. R8BHBC1( convolve13 )
  310. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  311. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  312. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  313. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  314. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  315. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  316. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  317. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  318. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]) +
  319. flt[ 9 ] * ( rp[ 19 ] + rp[ -19 ]) +
  320. flt[ 10 ] * ( rp[ 21 ] + rp[ -21 ]) +
  321. flt[ 11 ] * ( rp[ 23 ] + rp[ -23 ]) +
  322. flt[ 12 ] * ( rp[ 25 ] + rp[ -25 ]);
  323. R8BHBC2
  324. R8BHBC1( convolve14 )
  325. flt[ 0 ] * ( rp[ 1 ] + rp[ -1 ]) +
  326. flt[ 1 ] * ( rp[ 3 ] + rp[ -3 ]) +
  327. flt[ 2 ] * ( rp[ 5 ] + rp[ -5 ]) +
  328. flt[ 3 ] * ( rp[ 7 ] + rp[ -7 ]) +
  329. flt[ 4 ] * ( rp[ 9 ] + rp[ -9 ]) +
  330. flt[ 5 ] * ( rp[ 11 ] + rp[ -11 ]) +
  331. flt[ 6 ] * ( rp[ 13 ] + rp[ -13 ]) +
  332. flt[ 7 ] * ( rp[ 15 ] + rp[ -15 ]) +
  333. flt[ 8 ] * ( rp[ 17 ] + rp[ -17 ]) +
  334. flt[ 9 ] * ( rp[ 19 ] + rp[ -19 ]) +
  335. flt[ 10 ] * ( rp[ 21 ] + rp[ -21 ]) +
  336. flt[ 11 ] * ( rp[ 23 ] + rp[ -23 ]) +
  337. flt[ 12 ] * ( rp[ 25 ] + rp[ -25 ]) +
  338. flt[ 13 ] * ( rp[ 27 ] + rp[ -27 ]);
  339. R8BHBC2
  340. #undef R8BHBC1
  341. #undef R8BHBC2
  342. };
  343. // ---------------------------------------------------------------------------
  344. } // namespace r8b
  345. #endif // R8B_CDSPHBDOWNSAMPLER_INCLUDED