CDSPSincFilterGen.h

//$ nobt
//$ nocpp

/**
 * @file CDSPSincFilterGen.h
 *
 * @brief Sinc function-based FIR filter generator class.
 *
 * This file includes the CDSPSincFilterGen class implementation that
 * generates FIR filters.
 *
 * r8brain-free-src Copyright (c) 2013-2022 Aleksey Vaneev
 * See the "LICENSE" file for license.
 */

#ifndef R8B_CDSPSINCFILTERGEN_INCLUDED
#define R8B_CDSPSINCFILTERGEN_INCLUDED

#include "r8bbase.h"

namespace r8b {

/**
 * @brief Sinc function-based FIR filter generator class.
 *
 * Structure that holds state used to perform generation of sinc functions of
 * various types, windowed by the Blackman window by default (but the window
 * function can be changed if necessary).
 */

class CDSPSincFilterGen
{
public:
	double Len2; ///< Required half filter kernel's length in samples (can be
		///< a fractional value). Final physical kernel length will be
		///< provided in the KernelLen variable. Len2 should be >= 2.
	double Len2i; ///< = 1.0 / Len2, initialized and used by some window
		///< functions for optimization (should not be initialized by the
		///< caller).
	int KernelLen; ///< Resulting length of the filter kernel, this variable
		///< is set after the call to one of the "init" functions.
	int fl2; ///< Internal "half kernel length" value. This value can be used
		///< as filter's latency in samples (taps), this variable is set after
		///< the call to one of the "init" functions.

	union
	{
		struct
		{
			double Freq1; ///< Required corner circular frequency 1 [0; pi].
				///< Used only in the generateBand() function.
			double Freq2; ///< Required corner circular frequency 2 [0; pi].
				///< Used only in the generateBand() function. The range
				///< [Freq1; Freq2] defines a pass band for the generateBand()
				///< function.
		};

		struct
		{
			double FracDelay; ///< Fractional delay in the range [0; 1], used
				///< only in the generateFrac() function. Note that the
				///< FracDelay parameter is actually inversed. At 0.0 value it
				///< produces 1 sample delay (with the latency equal to fl2),
				///< at 1.0 value it produces 0 sample delay (with the latency
				///< equal to fl2 - 1).
		};
	};

	/**
	 * Window function type.
	 */

	enum EWindowFunctionType
	{
		wftCosine, ///< Generalized cosine window function. No parameters
			///< required. The "Power" parameter is optional.
		wftKaiser, ///< Kaiser window function. Requires the "Beta" parameter.
			///< The "Power" parameter is optional.
		wftGaussian ///< Gaussian window function. Requires the "Sigma"
			///< parameter. The "Power" parameter is optional.
	};

	typedef double( CDSPSincFilterGen :: *CWindowFunc )(); ///< Window
		///< calculation function pointer type.

	/**
	 * Function initializes *this structure for generation of a window
	 * function, odd-sized.
	 *
	 * @param WinType Window function type.
	 * @param Params Window function's parameters. If NULL, the table values
	 * may be used.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function. If "true", the power factor should be specified as the
	 * last value in the Params array. If Params is NULL, the table or default
	 * value of -1.0 (off) will be used.
	 */

	void initWindow( const EWindowFunctionType WinType = wftCosine,
		const double* const Params = NULL, const bool UsePower = false )
	{
		R8BASSERT( Len2 >= 2.0 );

		fl2 = (int) floor( Len2 );
		KernelLen = fl2 + fl2 + 1;

		setWindow( WinType, Params, UsePower, true );
	}

	/**
	 * Function initializes *this structure for generation of band-limited
	 * sinc filter kernel. The generateBand() function should be used to
	 * calculate the filter.
	 *
	 * @param WinType Window function type.
	 * @param Params Window function's parameters. If NULL, the table values
	 * may be used.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function. If "true", the power factor should be specified as the
	 * last value in the Params array. If Params is NULL, the table or default
	 * value of -1.0 (off) will be used.
	 */

	void initBand( const EWindowFunctionType WinType = wftCosine,
		const double* const Params = NULL, const bool UsePower = false )
	{
		R8BASSERT( Len2 >= 2.0 );

		fl2 = (int) floor( Len2 );
		KernelLen = fl2 + fl2 + 1;

		setWindow( WinType, Params, UsePower, true );
	}

	/**
	 * Function initializes *this structure for Hilbert transformation filter
	 * calculation. Freq1 and Freq2 variables are not used.
	 * The generateHilbert() function should be used to calculate the filter.
	 *
	 * @param WinType Window function type.
	 * @param Params Window function's parameters. If NULL, the table values
	 * may be used.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function. If "true", the power factor should be specified as the
	 * last value in the Params array. If Params is NULL, the table or default
	 * value of -1.0 (off) will be used.
	 */

	void initHilbert( const EWindowFunctionType WinType = wftCosine,
		const double* const Params = NULL, const bool UsePower = false )
	{
		R8BASSERT( Len2 >= 2.0 );

		fl2 = (int) floor( Len2 );
		KernelLen = fl2 + fl2 + 1;

		setWindow( WinType, Params, UsePower, true );
	}

	/**
	 * Function initializes *this structure for generation of full-bandwidth
	 * fractional delay sinc filter kernel. Freq1 and Freq2 variables are not
	 * used. The generateFrac() function should be used to calculate the
	 * filter.
	 *
	 * @param WinType Window function type.
	 * @param Params Window function's parameters. If NULL, the table values
	 * may be used.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function. If "true", the power factor should be specified as the
	 * last value in the Params array. If Params is NULL, the table or default
	 * value of -1.0 (off) will be used.
	 */

	void initFrac( const EWindowFunctionType WinType = wftCosine,
		const double* const Params = NULL, const bool UsePower = false )
	{
		R8BASSERT( Len2 >= 2.0 );

		fl2 = (int) ceil( Len2 );
		KernelLen = fl2 + fl2;

		setWindow( WinType, Params, UsePower, false, FracDelay );
	}

	/**
	 * @return The next "Hann" window function coefficient.
	 */

	double calcWindowHann()
	{
		return( 0.5 + 0.5 * w1.generate() );
	}

	/**
	 * @return The next "Hamming" window function coefficient.
	 */

	double calcWindowHamming()
	{
		return( 0.54 + 0.46 * w1.generate() );
	}

	/**
	 * @return The next "Blackman" window function coefficient.
	 */

	double calcWindowBlackman()
	{
		return( 0.42 + 0.5 * w1.generate() + 0.08 * w2.generate() );
	}

	/**
	 * @return The next "Nuttall" window function coefficient.
	 */

	double calcWindowNuttall()
	{
		return( 0.355768 + 0.487396 * w1.generate() +
			0.144232 * w2.generate() + 0.012604 * w3.generate() );
	}

	/**
	 * @return The next "Blackman-Nuttall" window function coefficient.
	 */

	double calcWindowBlackmanNuttall()
	{
		return( 0.3635819 + 0.4891775 * w1.generate() +
			0.1365995 * w2.generate() + 0.0106411 * w3.generate() );
	}

	/**
	 * @return The next "Kaiser" window function coefficient.
	 */

	double calcWindowKaiser()
	{
		const double n = 1.0 - sqr( wn * Len2i + KaiserLen2Frac );
		wn++;

		if( n <= 0.0 )
		{
			return( 0.0 );
		}

		return( besselI0( KaiserBeta * sqrt( n )) * KaiserMul );
	}

	/**
	 * @return The next "Gaussian" window function coefficient.
	 */

	double calcWindowGaussian()
	{
		const double f = exp( -0.5 * sqr( wn * GaussianSigmaI +
			GaussianSigmaFrac ));

		wn++;

		return( f );
	}

	/**
	 * Function calculates window function only.
	 *
	 * @param[out] op Output buffer, length = KernelLen.
	 * @param wfunc Window calculation function to use.
	 */

	void generateWindow( double* op,
		CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
	{
		op += fl2;
		double* op2 = op;

		int l = fl2;

		if( Power < 0.0 )
		{
			*op = ( *this.*wfunc )();

			while( l > 0 )
			{
				const double v = ( *this.*wfunc )();

				op++;
				op2--;
				*op = v;
				*op2 = v;
				l--;
			}
		}
		else
		{
			*op = pow_a(( *this.*wfunc )(), Power );

			while( l > 0 )
			{
				const double v = pow_a(( *this.*wfunc )(), Power );

				op++;
				op2--;
				*op = v;
				*op2 = v;
				l--;
			}
		}
	}

	/**
	 * Function calculates band-limited windowed sinc function-based filter
	 * kernel.
	 *
	 * @param[out] op Output buffer, length = KernelLen.
	 * @param wfunc Window calculation function to use.
	 */

	void generateBand( double* op,
		CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
	{
		CSineGen f2( Freq2, 0.0, 1.0 / R8B_PI );
		f2.generate();

		op += fl2;
		double* op2 = op;
		const double pw = Power;
		int t = 1;

		if( Freq1 < 2.3e-13 )
		{
			if( pw < 0.0 )
			{
				*op = Freq2 * ( *this.*wfunc )() / R8B_PI;

				while( t <= fl2 )
				{
					const double v = f2.generate() * ( *this.*wfunc )() / t;
					op++;
					op2--;
					*op = v;
					*op2 = v;
					t++;
				}
			}
			else
			{
				*op = Freq2 * pow_a(( *this.*wfunc )(), pw ) / R8B_PI;

				while( t <= fl2 )
				{
					const double v = f2.generate() *
						pow_a(( *this.*wfunc )(), pw ) / t;

					op++;
					op2--;
					*op = v;
					*op2 = v;
					t++;
				}
			}
		}
		else
		{
			CSineGen f1( Freq1, 0.0, 1.0 / R8B_PI );
			f1.generate();

			if( pw < 0.0 )
			{
				*op = ( Freq2 - Freq1 ) * ( *this.*wfunc )() / R8B_PI;

				while( t <= fl2 )
				{
					const double v = ( f2.generate() - f1.generate() ) *
						( *this.*wfunc )() / t;

					op++;
					op2--;
					*op = v;
					*op2 = v;
					t++;
				}
			}
			else
			{
				*op = ( Freq2 - Freq1 ) *
					pow_a(( *this.*wfunc )(), pw ) / R8B_PI;

				while( t <= fl2 )
				{
					const double v = ( f2.generate() - f1.generate() ) *
						pow_a(( *this.*wfunc )(), pw ) / t;

					op++;
					op2--;
					*op = v;
					*op2 = v;
					t++;
				}
			}
		}
	}

	/**
	 * Function calculates windowed Hilbert transformer filter kernel.
	 *
	 * @param[out] op Output buffer, length = KernelLen.
	 * @param wfunc Window calculation function to use.
	 */

	void generateHilbert( double* op,
		CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman )
	{
		static const double fvalues[ 2 ] = { 0.0, 2.0 / R8B_PI };
		op += fl2;
		double* op2 = op;

		( *this.*wfunc )();
		*op = 0.0;

		int t = 1;

		if( Power < 0.0 )
		{
			while( t <= fl2 )
			{
				const double v = fvalues[ t & 1 ] * ( *this.*wfunc )() / t;
				op++;
				op2--;
				*op = v;
				*op2 = -v;
				t++;
			}
		}
		else
		{
			while( t <= fl2 )
			{
				const double v = fvalues[ t & 1 ] *
					pow_a(( *this.*wfunc )(), Power ) / t;

				op++;
				op2--;
				*op = v;
				*op2 = -v;
				t++;
			}
		}
	}

	/**
	 * Function calculates windowed fractional delay filter kernel.
	 *
	 * @param[out] op Output buffer, length = KernelLen.
	 * @param wfunc Window calculation function to use.
	 * @param opinc Output buffer increment, in "op" elements.
	 */

	void generateFrac( double* op,
		CWindowFunc wfunc = &CDSPSincFilterGen :: calcWindowBlackman,
		const int opinc = 1 )
	{
		R8BASSERT( opinc != 0 );

		const double pw = Power;
		const double fd = FracDelay;
		int t = -fl2;

		if( t + fd < -Len2 )
		{
			( *this.*wfunc )();
			*op = 0.0;
			op += opinc;
			t++;
		}

		double f = sin( fd * R8B_PI ) / R8B_PI;

		if(( t & 1 ) != 0 )
		{
			f = -f;
		}

		int IsZeroX = ( fabs( fd - 1.0 ) < 2.3e-13 );
		int mt = 0 - IsZeroX;
		IsZeroX = ( IsZeroX || fabs( fd ) < 2.3e-13 );

		if( pw < 0.0 )
		{
			while( t < mt )
			{
				*op = f * ( *this.*wfunc )() / ( t + fd );
				op += opinc;
				t++;
				f = -f;
			}

			if( IsZeroX ) // t+FracDelay==0
			{
				*op = ( *this.*wfunc )();
			}
			else
			{
				*op = f * ( *this.*wfunc )() / fd; // t==0
			}

			mt = fl2 - 2;

			while( t < mt )
			{
				op += opinc;
				t++;
				f = -f;
				*op = f * ( *this.*wfunc )() / ( t + fd );
			}

			op += opinc;
			t++;
			f = -f;
			const double ut = t + fd;
			*op = ( ut > Len2 ? 0.0 : f * ( *this.*wfunc )() / ut );
		}
		else
		{
			while( t < mt )
			{
				*op = f * pow_a(( *this.*wfunc )(), pw ) / ( t + fd );
				op += opinc;
				t++;
				f = -f;
			}

			if( IsZeroX ) // t+FracDelay==0
			{
				*op = pow_a(( *this.*wfunc )(), pw );
			}
			else
			{
				*op = f * pow_a(( *this.*wfunc )(), pw ) / fd; // t==0
			}

			mt = fl2 - 2;

			while( t < mt )
			{
				op += opinc;
				t++;
				f = -f;
				*op = f * pow_a(( *this.*wfunc )(), pw ) / ( t + fd );
			}

			op += opinc;
			t++;
			f = -f;
			const double ut = t + FracDelay;
			*op = ( ut > Len2 ? 0.0 :
				f * pow_a(( *this.*wfunc )(), pw ) / ut );
		}
	}

private:
	double Power; ///< The power factor used to raise the window function.
		///< Equals a negative value if the power factor should not be used.
	CSineGen w1; ///< Cosine wave 1 for window function.
	CSineGen w2; ///< Cosine wave 2 for window function.
	CSineGen w3; ///< Cosine wave 3 for window function.

	union
	{
		struct
		{
			double KaiserBeta; ///< Kaiser window function's "Beta"
				///< coefficient.
			double KaiserMul; ///< Kaiser window function's divisor, inverse.
			double KaiserLen2Frac; ///< Equals FracDelay / Len2.
		};

		struct
		{
			double GaussianSigmaI; ///< Gaussian window function's "Sigma"
				///< coefficient, inverse.
			double GaussianSigmaFrac; ///< Equals FracDelay / GaussianSigma.
		};
	};

	int wn; ///< Window function integer position. 0 - center of the window
		///< function. This variable may not be used by some window functions.

	/**
	 * Function initializes Kaiser window function calculation. The FracDelay
	 * variable should be initialized when using this window function.
	 *
	 * @param Params Function parameters. If NULL, the default values will be
	 * used. If not NULL, the first parameter should specify the "Beta" value.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function.
	 * @param IsCentered "True" if centered window should be used. This
	 * parameter usually equals to "false" for fractional delay filters only.
	 */

	void setWindowKaiser( const double* Params, const bool UsePower,
		const bool IsCentered )
	{
		wn = ( IsCentered ? 0 : -fl2 );

		if( Params == NULL )
		{
			KaiserBeta = 9.5945013206755156;
			Power = ( UsePower ? 1.9718457932433306 : -1.0 );
		}
		else
		{
			KaiserBeta = clampr( Params[ 0 ], 1.0, 350.0 );
			Power = ( UsePower ? fabs( Params[ 1 ]) : -1.0 );
		}

		KaiserMul = 1.0 / besselI0( KaiserBeta );
		Len2i = 1.0 / Len2;
		KaiserLen2Frac = FracDelay * Len2i;
	}

	/**
	 * Function initializes Gaussian window function calculation. The FracDelay
	 * variable should be initialized when using this window function.
	 *
	 * @param Params Function parameters. If NULL, the table values will be
	 * used. If not NULL, the first parameter should specify the "Sigma"
	 * value.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function.
	 * @param IsCentered "True" if centered window should be used. This
	 * parameter usually equals to "false" for fractional delay filters only.
	 */

	void setWindowGaussian( const double* Params, const bool UsePower,
		const bool IsCentered )
	{
		wn = ( IsCentered ? 0 : -fl2 );

		if( Params == NULL )
		{
			GaussianSigmaI = 1.0;
			Power = -1.0;
		}
		else
		{
			GaussianSigmaI = clampr( fabs( Params[ 0 ]), 1e-1, 100.0 );
			Power = ( UsePower ? fabs( Params[ 1 ]) : -1.0 );
		}

		GaussianSigmaI *= Len2;
		GaussianSigmaI = 1.0 / GaussianSigmaI;
		GaussianSigmaFrac = FracDelay * GaussianSigmaI;
	}

	/**
	 * Function initializes calculation of window function of the specified
	 * type.
	 *
	 * @param WinType Window function type.
	 * @param Params Window function's parameters. If NULL, the table values
	 * may be used.
	 * @param UsePower "True" if the power factor should be used to raise the
	 * window function. If "true", the power factor should be specified as the
	 * last value in the Params array. If Params is NULL, the table or default
	 * value of -1.0 (off) will be used.
	 * @param IsCentered "True" if centered window should be used. This
	 * parameter usually equals to "false" for fractional delay filters only.
	 * @param UseFracDelay Fractional delay to use.
	 */

	void setWindow( const EWindowFunctionType WinType,
		const double* const Params, const bool UsePower,
		const bool IsCentered, const double UseFracDelay = 0.0 )
	{
		FracDelay = UseFracDelay;

		if( WinType == wftCosine )
		{
			if( IsCentered )
			{
				w1.init( R8B_PI / Len2, R8B_PId2 );
				w2.init( R8B_2PI / Len2, R8B_PId2 );
				w3.init( R8B_3PI / Len2, R8B_PId2 );
			}
			else
			{
				const double step1 = R8B_PI / Len2;
				w1.init( step1, R8B_PId2 - step1 * fl2 + step1 * FracDelay );

				const double step2 = R8B_2PI / Len2;
				w2.init( step2, R8B_PId2 - step2 * fl2 + step2 * FracDelay );

				const double step3 = R8B_3PI / Len2;
				w3.init( step3, R8B_PId2 - step3 * fl2 + step3 * FracDelay );
			}

			Power = ( UsePower && Params != NULL ? Params[ 0 ] : -1.0 );
		}
		else
		if( WinType == wftKaiser )
		{
			setWindowKaiser( Params, UsePower, IsCentered );
		}
		else
		if( WinType == wftGaussian )
		{
			setWindowGaussian( Params, UsePower, IsCentered );
		}
	}
};

} // namespace r8b

#endif // R8B_CDSPSINCFILTERGEN_INCLUDED