robert@464: /*
robert@464:  *
robert@464:  * Simple 1-Dimensional Waveguide
robert@464:  *
robert@464:  */
robert@464: 
robert@464: #include "Waveguide.h"
robert@464: #include "../include/Utilities.h"
robert@464: #include <rtdk.h>
robert@464: #include <cmath>
robert@464: #include <stdio.h>
robert@464: #include <cstdlib>
robert@464: 
robert@464: #define DECAY 0.995//0.999
robert@464: #define DAMPING 0.01//0.05
robert@464: 
robert@464: // TODO: make damping and decay parametrisable
robert@464: 
robert@464: Waveguide::Waveguide()	{
robert@464: 
robert@464: 	// initialize variables
robert@464: 	a1_lp = 0;
robert@464: 	a2_lp = 0;
robert@464: 	b0_lp = 0;
robert@464: 	b1_lp = 0;
robert@464: 	b2_lp = 0;
robert@464: 	_dt = 1.0/44100.0;
robert@464: 	setFrequency(440);
robert@464: 	updateFilterCoeffs(8000);
robert@464: 	_filterReadPtr=0;
robert@464: 	for(int i=0;i<FILTER_BUFFER_SIZE;i++)	{
robert@464: 		_filterBuffer_x[i] = 0;
robert@464: 		_filterBuffer_y[i] = 0;
robert@464: 	}
robert@464: 	for(int i=0;i<WG_BUFFER_SIZE;i++)	{
robert@464: 		_buffer[i] = 0;
robert@464: 	}
robert@464: 	_lastX = 0;
robert@464: 	_lastY = 0;
robert@464: 	_readPtr = 0;
robert@464: 
robert@464: }
robert@464: 
robert@464: void Waveguide::setup()	{
robert@464: 
robert@464: }
robert@464: 
robert@464: float Waveguide::update(float in)	{
robert@464: 
robert@464: 	// 1. advance delay buffer read pointer
robert@464: 
robert@464: 	if(++_readPtr>=WG_BUFFER_SIZE)
robert@464: 		_readPtr=0;
robert@464: 
robert@464: 	// 2. write input into buffer
robert@464: 
robert@464: 	_buffer[_readPtr] = in;
robert@464: 
robert@464: 	// 3. read delayed sample from buffer
robert@464: 
robert@464: 	float out = _buffer[(_readPtr-_periodInSamples+WG_BUFFER_SIZE)%WG_BUFFER_SIZE];
robert@464: 
robert@464: 	// 4. apply damping (low-pass) filter to output
robert@464: 
robert@464: 	if(++_filterReadPtr>=FILTER_BUFFER_SIZE)
robert@464: 			_filterReadPtr=0;
robert@464: 
robert@464: 	out = b0_lp*out +
robert@464: 					b1_lp*_filterBuffer_x[(_filterReadPtr-1+FILTER_BUFFER_SIZE)%FILTER_BUFFER_SIZE] +
robert@464: 					b2_lp*_filterBuffer_x[(_filterReadPtr-2+FILTER_BUFFER_SIZE)%FILTER_BUFFER_SIZE] -
robert@464: 					a1_lp*_filterBuffer_y[(_filterReadPtr-1+FILTER_BUFFER_SIZE)%FILTER_BUFFER_SIZE] -
robert@464: 					a2_lp*_filterBuffer_y[(_filterReadPtr-2+FILTER_BUFFER_SIZE)%FILTER_BUFFER_SIZE];
robert@464: 
robert@464: 	// 5. Simple high-pass filter to block DC-offset
robert@464: 	// y[n] = x[n] - x[n-1] + a * y[n-1]
robert@464: 	float gain = 0.9999;
robert@464: 	float temp = out;
robert@464: 	out = out - _lastX + gain * _lastY;
robert@464: 	_lastY = out;
robert@464: 	_lastX = temp;
robert@464: 
robert@464: 	// 6. Apply intensity damping
robert@464: 	out *= DECAY;
robert@464: 
robert@464: 	_filterBuffer_x[_filterReadPtr] = in;
robert@464: 	_filterBuffer_y[_filterReadPtr] = out;
robert@464: 
robert@464: 	return out;
robert@464: 
robert@464: }
robert@464: 
robert@464: void Waveguide::setFrequency(float frequency)	{
robert@464: 
robert@464: 	// NB: currently no interpolation, so may not be ideal for dynamically changing waveguide frequency
robert@464: 	_periodInMilliseconds = 1000.0/frequency;
robert@464: 	_periodInSamples = (int)(_periodInMilliseconds * 44.1);
robert@464: 
robert@464: }
robert@464: 
robert@464: void Waveguide::updateFilterCoeffs(float frequency)	{
robert@464: 
robert@464: 	// FIXME: Butterworth filter doesn't work very well,
robert@464: 	//		  using simple FIR in the meantime
robert@464: 
robert@464: 	a1_lp = 0;
robert@464: 	a2_lp = 0;
robert@464: 	b0_lp = 1.0 - DAMPING;
robert@464: 	b1_lp = DAMPING;
robert@464: 	b2_lp = 0;
robert@464: 
robert@464: 	/*
robert@464: 	// 'w' for sake of resembling lower-case 'omega'
robert@464: 	float w = 2.0 * M_PI * frequency;
robert@464: 	float t = _dt;
robert@464: 	// The Q for a 2nd-order Butterworth is sqrt(2)/2
robert@464: 	float q = 0.707;//sqrt(2.0)/2.0;
robert@464: 
robert@464: 	// low-pass filter coefficients
robert@464: 	float a0_lp = w*w*t*t + 2*(w/q)*t + 4.0;
robert@464: 	float k = 1.0/a0_lp;
robert@464: 	a1_lp = (2.0*w*w*t*t - 8.0) * k;
robert@464: 	a2_lp = (4.0 - (w/q)*2.0*t + w*w*t*t) * k;
robert@464: 	b0_lp = (w*w*t*t) * k;
robert@464: 	b1_lp = (2.0*w*w*t*t) * k;
robert@464: 	b2_lp = (w*w*t*t) * k;
robert@464: 	*/
robert@464: 
robert@464: }