snd_dsp.cpp

hl2 source code. Do not use it illegal.

// snd_dsp.c -- audio processing routines

#define WIN32_LEAN_AND_MEAN
#pragma warning(push, 1)
#include <windows.h>
#include <mmsystem.h>
#pragma warning(pop)

#include "tier0/dbg.h"
#include "sound.h"
#include "sound_private.h"
#include "soundflags.h"
#include "snd_device.h"
#include "measure_section.h"
#include "iprediction.h"
#include "snd_mix_buf.h"
#include "snd_env_fx.h"
#include "snd_channels.h"
#include "snd_audio_source.h"
#include "snd_convars.h"
#include "SoundService.h"
#include "commonmacros.h"

#include "mathlib.h"


// hard clip input value to -32767 <= y <= 32767

// #define CLIP(x) ((x) > 32767 ? 32767 : ((x) < -32767 ? -32767 : (x)))

#define SIGN(d)				((d)<0?-1:1)

#define ABS(a)	abs(a)

#define MSEC_TO_SAMPS(a)	(((a)*SOUND_DMA_SPEED) / 1000)		// convert milliseconds to # samples in equivalent time
#define SEC_TO_SAMPS(a)		((a)*SOUND_DMA_SPEED)				// conver seconds to # samples in equivalent time

//#define CLIP_DSP(x) ((x) > 32767 ? 32767 : ((x) < -32767 ? -32767 : (x)))

#define CLIP_DSP(x) (x)

extern bool SURROUND_ON;

#define SOUND_MS_PER_FT	1			// sound travels approx 1 foot per millisecond
#define ROOM_MAX_SIZE	1000		// max size in feet of room simulation for dsp

//===============================================================================
//
// Digital Signal Processing algorithms for audio FX.
//
// KellyB 2/18/03
//===============================================================================

// Performance notes:

// DSP processing should take no more than 3ms total time per frame to remain on par with hl1
// Assume a min frame rate of 24fps = 42ms per frame
// at 24fps, to maintain 44.1khz output rate, we must process about 1840 mono samples per frame.
// So we must process 1840 samples in 3ms.

// on a 1Ghz CPU (mid-low end CPU) 3ms provides roughly 3,000,000 cycles.
// Thus we have 3e6 / 1840 = 1630 cycles per sample.  

#define PBITS			12					// parameter bits
#define PMAX			((1 << PBITS)-1)	// parameter max size

// crossfade from y2 to y1 at point r (0 < r < PMAX)

#define XFADE(y1,y2,r)	(((y1) * (r)) >> PBITS) + (((y2) * (PMAX - (r))) >> PBITS);

#define XFADEF(y1,y2,r)	(((y1) * (r)) / (float)(PMAX)) + (((y2) * (PMAX - (r))) / (float)(PMAX));

/////////////////////
// dsp helpers
/////////////////////

// dot two integer vectors of length M+1
// M is filter order, h is filter vector, w is filter state vector

inline int dot ( int M, int *h, int *w )
{
	int i;
	int y;

	for (y = 0, i = 0; i <= M; i++)
		y += ( h[i] * w[i] ) >> PBITS;
	return y;
}

// delay array w[] by D samples
// w[0] = input, w[D] = output
// practical for filters, but not for large values of D

inline void delay ( int D, int *w )
{
	int i;
	for (i = D; i >= 1; i--)		// reverse order updating
		w[i] = w[i-1];
}

// circular wrap of pointer p, relative to array w
// D delay line size in samples w[0...D]
// w delay line buffer pointer, dimension D+1
// p circular pointer

inline void wrap (int D, int *w, int **p)
{
	if ( *p > w + D )
		*p -= D + 1;		// when *p = w + D + 1, it wraps around to *p = w
	
	if ( *p < w )
		*p += D + 1;		// when *p = w - 1, it wraps around to *p = w + D
}

// simple averaging filter for performance - a[] is 0, b[] is 1, L is # of samples to average

inline int avg_filter ( int  M, int *a, int L, int *b, int *w, int x )
{
	int i;
	int y = 0;
	
	
	w[0] = x;

	// output adder

	switch (L)
	{
	default:
	case 12: y += w[12];
	case 11: y += w[11];
	case 10: y += w[10];
	case 9:  y += w[9];
	case 8:  y += w[8];
	case 7:  y += w[7];
	case 6:  y += w[6];
	case 5:  y += w[5];
	case 4:  y += w[4];
	case 3:  y += w[3];
	case 2:  y += w[2];
	case 1:  y += w[1];
	case 0:  y += w[0];
	}
	
	for (i = L; i >= 1; i--)		// reverse update internal state
		w[i] = w[i-1];

	switch (L)
	{
	default:
	case 12: return y / 13;
	case 11: return y / 12;
	case 10: return y / 11;
	case 9:  return y / 10;
	case 8:  return y / 9;
	case 7:  return y >> 3;
	case 6:  return y / 7;
	case 5:  return y / 6;
	case 4:  return y / 5;
	case 3:  return y >> 2;
	case 2:  return y / 3;
	case 1:  return y >> 1;
	case 0:  return y;
	}

	return y;						// current output sample
}

// IIR filter, cannonical form
//  returns single sample y for current input value x
//  x is input sample
//	w = internal state vector, dimension max(M,L) + 1
//  L, M numerator and denominator filter orders
//  a,b are M+1 dimensional arrays of filter params
// 
//  for M = 4:
//
//                1     w0(n)   b0	 
//  x(n)--->(+)--(*)-----.------(*)->(+)---> y(n)
//           ^         |            ^
//           |     [Delay d]        |
//           |         |            |
//           |    -a1  |W1     b1   |
//			 ----(*)---.------(*)----
//           ^         |            ^
//           |     [Delay d]        |
//           |         |            |
//           |    -a2  |W2     b2   |
//			 ----(*)---.------(*)----
//           ^         |            ^
//           |     [Delay d]        |
//           |         |            |
//           |    -a3  |W3     b3   |
//			 ----(*)---.------(*)----
//           ^         |            ^
//           |     [Delay d]        |
//           |         |            |
//           |    -a4  |W4     b4   |
//			 ----(*)---.------(*)----
//
//	for each input sample x, do:
//			w0 = x - a1*w1 - a2*w2 - ... aMwM
//			y = b0*w0 + b1*w1 + ...bL*wL
//			wi = wi-1, i = K, K-1, ..., 1


inline int iir_filter ( int  M, int *a, int L, int *b, int *w, int x )
{
	int K, i;
	int y;
	int x0;					
	
	if (M == 0)
		return (avg_filter( M, a, L, b, w, x));
	
	y = 0;
	x0 = x;

	K = max ( M, L );				
	
	// for (i = 1; i <= M; i++)		// input adder
	//	w[0] -= ( a[i] * w[i] ) >> PBITS;

	// M is clamped between 1 and FLT_M
	// change this switch statement if FLT_M changes!

	switch (M)
	{
	case 12: x0 -= ( a[12] * w[12] ) >> PBITS;
	case 11: x0 -= ( a[11] * w[11] ) >> PBITS;
	case 10: x0 -= ( a[10] * w[10] ) >> PBITS;
	case 9:  x0 -= ( a[9] * w[9] ) >> PBITS;
	case 8:  x0 -= ( a[8] * w[8] ) >> PBITS;
	case 7:  x0 -= ( a[7] * w[7] ) >> PBITS;
	case 6:  x0 -= ( a[6] * w[6] ) >> PBITS;
	case 5:  x0 -= ( a[5] * w[5] ) >> PBITS;
	case 4:  x0 -= ( a[4] * w[4] ) >> PBITS;
	case 3:  x0 -= ( a[3] * w[3] ) >> PBITS;
	case 2:  x0 -= ( a[2] * w[2] ) >> PBITS;
	default:
	case 1:  x0 -= ( a[1] * w[1] ) >> PBITS;
	}

	w[0] = x0;

	//for (i = 0; i <= L; i++)		// output adder
	//	y += ( b[i] * w[i] ) >> PBITS;

	switch (L)
	{
	case 12: y += ( b[12] * w[12] ) >> PBITS;
	case 11: y += ( b[11] * w[11] ) >> PBITS;
	case 10: y += ( b[10] * w[10] ) >> PBITS;
	case 9:  y += ( b[9] * w[9] ) >> PBITS;
	case 8:  y += ( b[8] * w[8] ) >> PBITS;
	case 7:  y += ( b[7] * w[7] ) >> PBITS;
	case 6:  y += ( b[6] * w[6] ) >> PBITS;
	case 5:  y += ( b[5] * w[5] ) >> PBITS;
	case 4:  y += ( b[4] * w[4] ) >> PBITS;
	case 3:  y += ( b[3] * w[3] ) >> PBITS;
	case 2:  y += ( b[2] * w[2] ) >> PBITS;
	default:
	case 1:  y += ( b[1] * w[1] ) >> PBITS;
	case 0:  y += ( b[0] * w[0] ) >> PBITS;
	}
	
	for (i = K; i >= 1; i--)		// reverse update internal state
		w[i] = w[i-1];

	return y;						// current output sample
}

// IIR filter, cannonical form, using dot product and delay implementation
// (may be easier to optimize this routine.)

inline int iir_filter2 ( int  M, int *a, int L, int *b, int *w, int x )
{
	int K;
	int y;

	K = max ( M, L );				// K = max (M, L)
	w[0] = 0;						// needed for dot (M, a, w)
	
	w[0] = x - dot ( M, a, w );		// input adder
	y = dot ( L, b, w );			// output adder
	
	delay ( K, w );					// update delay line
	
	return y;						// current output sample
}


// fir filter - no feedback = high stability but also may be more expensive computationally

inline int fir_filter ( int M, int *h, int *w, int x )
{
	int i;
	int y;

	w[0] = x;

	for ( y = 0, i = 0; i <= M; i++ )
		y += h[i] * w[i];

	for ( i = M; i >= -1; i-- )
		w[i] = w[i-1];

	return y;
}

// fir filter, using dot product and delay implementation

inline int fir_filter2 ( int M, int *h, int *w, int x )
{
	int y;

	w[0] = x;

	y = dot ( M, h, w );

	delay ( M, w );

	return y;
}


// tap - i-th tap of circular delay line buffer
// D delay line size in samples
// w delay line buffer pointer, of dimension D+1
// p circular pointer
// t = 0...D

int tap ( int D, int *w, int *p, int t )
{
	return w[(p - w + t) % (D + 1)];
}

// tapi - interpolated tap output of a delay line
//        interpolates sample between adjacent samples in delay line for 'frac' part of delay
// D delay line size in samples
// w delay line buffer pointer, of dimension D+1
// p circular pointer
// t - delay tap integer value 0...D.  (complete delay is t.frac )
// frac - varying 16 bit fractional delay value 0...32767 (normalized to 0.0 - 1.0)

inline int tapi ( int D, int *w, int *p, int t, int frac )
{
	int i, j;
	int si, sj;

	i = t;					// tap value, interpolate between adjacent samples si and sj
	j = (i + 1) % (D+1);	// if i = D, then j = 0; otherwise, j = i + 1
	
	si = tap( D, w, p, i );	// si(n) = x(n - i)
	sj = tap( D, w, p, j );	// sj(n) = x(n - j)

	return si + (((frac) * (sj - si) ) >> 16);
}

// circular delay line, D-fold delay
// D delay line size in samples w[0..D]
// w delay line buffer pointer, dimension D+1
// p circular pointer

inline void cdelay ( int D, int *w, int **p )
{
	(*p)--;					// decrement pointer and wrap modulo (D+1)
	wrap ( D, w, p );		// when *p = w-1, it wraps around to *p = w+D
}

// plain reverberator with circular delay line
// D delay line size in samples
// t tap from this location - <= D
// w delay line buffer pointer of dimension D+1
// p circular pointer, must be init to &w[0] before first call
// a feedback value, 0-PMAX (normalized to 0.0-1.0)
// b gain
// x input sample

//                    w0(n)   b	 
//  x(n)--->(+)--------.-----(*)-> y(n)
//           ^         |            
//           |     [Delay d]        
//           |         |            
//           |    a    |Wd(n)       
//			 ----(*)---.

inline int dly_plain ( int D, int t, int *w, int **p, int a, int b, int x )
{
	int y, sD;

	sD = tap ( D, w, *p, t );		// Tth tap delay output
	y = x + (( a * sD ) >> PBITS);		// filter output
	**p = y;						// delay input
	cdelay ( D, w, p );				// update delay line

	return ( (y * b) >> PBITS );
}

// straight delay line
//
// D delay line size in samples
// t tap from this location - <= D
// w delay line buffer pointer of dimension D+1
// p circular pointer, must be init to &w[0] before first call
// x input sample
//                    
//  x(n)--->[Delay d]---> y(n)
//
 
inline int dly_linear ( int D, int t, int *w, int **p, int x )
{
	int y;

	y = tap ( D, w, *p, t );		// Tth tap delay output
	**p = x;						// delay input
	cdelay ( D, w, p );				// update delay line

	return ( y );
}

// lowpass reverberator, replace feedback multiplier 'a' in 
// plain reverberator with a low pass filter
// D delay line size in samples
// t tap from this location - <= D
// w delay line buffer pointer of dimension D+1
// p circular pointer, must be init to &w[0] before first call
// a feedback gain
// b output gain
// M filter order
// bf filter numerator, 0-PMAX (normalized to 0.0-1.0), M+1 dimensional
// af filter denominator, 0-PMAX (normalized to 0.0-1.0), M+1 dimensional
// vf filter state, M+1 dimensional
// x input sample
//                    w0(n)    	   b
//  x(n)--->(+)--------------.----(*)--> y(n)
//           ^               |            
//           |           [Delay d]        
//           |               |            
//           |  a            |Wd(n)       
//			 --(*)--[Filter])-

int dly_lowpass ( int D, int t, int *w, int **p, int a, int b, int M, int *af, int L, int *bf, int *vf, int x )
{
	int y, sD;

	sD = tap ( D, w, *p, t );										// delay output is filter input
	y = x + ((iir_filter ( M, af, L, bf, vf, sD ) * a) >> PBITS);	// filter output with gain
	**p = y;														// delay input
	cdelay ( D, w, p );												// update delay line

	return ( (y * b) >> PBITS ); // output with gain
}

// allpass reverberator with circular delay line
// D delay line size in samples
// t tap from this location - <= D
// w delay line buffer pointer of dimension D+1
// p circular pointer, must be init to &w[0] before first call
// a feedback value, 0-PMAX (normalized to 0.0-1.0)
// b gain

//                    w0(n)   -a	     b
//  x(n)--->(+)--------.-----(*)-->(+)--(*)-> y(n)
//           ^         |            ^
//           |     [Delay d]        |
//           |         |            |
//           |    a    |Wd(n)       |
//			 ----(*)---.-------------
//
//	for each input sample x, do:
//		w0 = x + a*Wd
//		y = -a*w0 + Wd
//		delay (d, W) - w is the delay buffer array
//
// or, using circular delay, for each input sample x do:
//
//		Sd = tap (D,w,p,D)
//		S0 = x + a*Sd
//		y = -a*S0 + Sd
//		*p = S0
//		cdelay(D, w, &p)		

inline int dly_allpass ( int D, int t, int *w, int **p, int a, int b, int x ) 
{
	int y, s0, sD;
	
	sD = tap ( D, w, *p, t );			// Dth tap delay output
	s0 = x + (( a * sD ) >> PBITS);

	y = ( ( -a * s0 ) >> PBITS ) + sD;	// filter output
	**p = s0;							// delay input
	cdelay ( D, w, p );					// update delay line

	return ( (y * b) >> PBITS );
}


///////////////////////////////////////////////////////////////////////////////////
// fixed point math for real-time wave table traversing, pitch shifting, resampling
///////////////////////////////////////////////////////////////////////////////////

#define FIX20_BITS			20									// 20 bits of fractional part
#define FIX20_SCALE			(1 << FIX20_BITS)

#define FIX20_INTMAX		((1 << (32 - FIX20_BITS))-1)		// maximum step integer

#define FLOAT_TO_FIX20(a)	((int)((a) * (float)FIX20_SCALE))		// convert float to fixed point
#define INT_TO_FIX20(a)		(((int)(a)) << FIX20_BITS)			// convert int to fixed point
#define FIX20_TO_FLOAT(a)	((float)(a) / (float)FIX20_SCALE)	// convert fix20 to float