Commit 2e236b0d authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files

Moving PD and JetEngine into own class / files

parent 5e6819a9
......@@ -11,6 +11,7 @@ vista_use_package( ITAFFT REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAConvolution REQUIRED FIND_DEPENDENCIES ) # required for FIR filtering (uses block convolver)
vista_use_package( TBB REQUIRED FIND_DEPENDENCIES )
vista_use_package( SPLINE REQUIRED FIND_DEPENDENCIES )
vista_use_package( DspFilters REQUIRED FIND_DEPENDENCIES )
vista_find_package( libjson OPTIONAL QUIET )
vista_find_package( MKL OPTIONAL QUIET )
......@@ -28,6 +29,7 @@ include_directories( "include" )
# sources
set( ITADSPHeader
"include/ITADSP/PD/JetEngine.h"
"include/ITADSP/ThirdOctaveFilterbank/IIRBurg.h"
"include/ITADSPDefinitions.h"
"include/ITASIMOVariableDelayLine.h"
......@@ -43,6 +45,7 @@ set( ITADSPHeader
"include/ITAThirdOctaveFIRFilterGenerator.h"
)
set( ITADSPSources
"src/ITADSP/PD/JetEngine.cpp"
"src/ITADSP/ThirdOctaveFilterbank/IIRBurg.cpp"
"src/ITASIMOVariableDelayLine.cpp"
"src/ITASIMOVariableDelayLineBase.cpp"
......
/*
* ----------------------------------------------------------------
*
* ITA core libs
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_DSP_PD_JET_ENGINE
#define IW_ITA_DSP_PD_JET_ENGINE
#include <ITADSPDefinitions.h>
#include <ITAThirdOctaveFilterbank.h>
#include <VistaTools/VistaRandomNumberGenerator.h>
#include <DspFilters/Butterworth.h>
#include <tbb/concurrent_queue.h>
#include <atomic>
#include <algorithm>
#include <atomic>
#include <vector>
class CITAIIRFilterEngine;
class CITAThirdOctaveFIRFilterGenerator;
const int g_iIIRFilterOrder = 1;
const int g_iChannels = 1;
namespace ITADSP
{
namespace PD
{
//! Pure data patch implementation of a jet engine
/**
* Manual implementation of the "jet engine" pure data patch (combining "forced flame" and "turbine") from Andy Farnell's Designing Sound
*/
class ITA_DSP_API CJetEngine
{
public:
inline CJetEngine( double dSampleRate = 44100.0f, float fRPMInit = 1000.0f, bool bColdStart = true )
: m_dSampleRate( dSampleRate )
{
m_vfTurbineBaseModeFrequencies = { 3097.0f, 4495.0f, 5588.0f, 7414.0f, 11000.0f };
m_vfTurbineBaseModeAmplitudes = { 0.25f, 0.25f, 1.0f, 0.4f, 0.4f };
m_vfRPMRange = { 1000.0f, 5000.0f };
const float fCenterFrequency1 = 8000.0f;
const float fQ1 = 0.5f;
m_oForcedFlameBP1.setup( g_iIIRFilterOrder, float( dSampleRate ), fCenterFrequency1, fCenterFrequency1 / fQ1 );
const float fCutoffFrequency1 = 120.0f;
m_oForcedFlameHP1.setup( g_iIIRFilterOrder, float( dSampleRate ), fCutoffFrequency1 );
const float fLPCutoffFrequency1 = 11000.0f;
m_oJetEngineLP1.setup( g_iIIRFilterOrder, float( dSampleRate ), fLPCutoffFrequency1 );
const float fLPCutoffFrequency2 = 0.1f;
m_oJetEngineLP2.setup( 1, float( dSampleRate ), fLPCutoffFrequency2 );
oCtxAudio.vfTurbineModePhaseShift = { 0, 0, 0, 0, 0 };
oCtxAudio.vfTurbineModeFrequencies = { 0, 0, 0, 0, 0 };
SetRPM( fRPMInit );
if( !bColdStart )
{
// Run input low-pass for 3 secs (until 0.2 Hz LP died away on unit step function) to akkumulate
// of the default input, otherwise there is a run-up from zero to target RPM value
const float fNormalizedRPMInit = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMUpsample;
float* pfNormalizedRPMUpsampleAlias = &fNormalizedRPMUpsample;
for( int n = 0; n < ( int ) dSampleRate * 3; n++ )
{
fNormalizedRPMUpsample = fNormalizedRPMInit;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMUpsampleAlias );
}
}
};
inline virtual ~CJetEngine()
{
};
virtual void SetRPM( float fNewRPM );
virtual void Process( float* pfOutputBuffer, int iNumSamples );
protected:
// We'll protect data exchange between user and audio thread
// with atomics, but a request, swap & acknowledge approach with
// mutex / locking would be more elegant
struct CAudioContext
{
std::atomic< float > vfRPM;
std::vector< float > vfTurbineModeFrequencies;
std::vector< float > vfTurbineModePhaseShift;
} oCtxAudio;
float GetRPMNormalized( float fRPM ) const;
private:
VistaRandomNumberGenerator oRNG;
Dsp::SimpleFilter< Dsp::Butterworth::BandPass< g_iIIRFilterOrder >, g_iChannels > m_oForcedFlameBP1, m_oForcedFlameBP2, m_oForcedFlameBP3;
Dsp::SimpleFilter< Dsp::Butterworth::HighPass< g_iIIRFilterOrder >, g_iChannels > m_oForcedFlameHP1;
Dsp::SimpleFilter< Dsp::Butterworth::LowPass< 1 >, g_iChannels > m_oJetEngineLP1, m_oJetEngineLP2;
std::vector< float > m_vfTurbineBaseModeFrequencies, m_vfTurbineBaseModeAmplitudes, m_vfRPMRange;
float m_fTempSample;
double m_dSampleRate;
};
}
}
#endif // IW_ITA_DSP_PD_JET_ENGINE
#include <ITADSP/PD/JetEngine.h>
void ITADSP::PD::CJetEngine::SetRPM( float fRPM )
{
oCtxAudio.vfRPM = fRPM; // Direct write-through, no lock & swap implemented here
}
float ITADSP::PD::CJetEngine::GetRPMNormalized( float fRPM ) const
{
// Validate input range
assert( m_vfRPMRange.size() == 2 );
const float fValidRPM = std::min( std::max( m_vfRPMRange[ 0 ], fRPM ), m_vfRPMRange[ 1 ] );
const float fNormalizedRPM = std::max( 0.1f, fValidRPM / m_vfRPMRange[ 1 ] );
return fNormalizedRPM;
}
void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
{
/* Validate, normalize and up-sampled RPM and apply low-pass filter */
// Processed DSP coefficient with normalized input
const float fNormalizedRPMControl = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMControlUpscaled;
// We need this ref pointer for Dsp functions, but its ugly we dont use it apart from that
float* pfTempSampleAlias = &m_fTempSample;
float* pfNormalizedRPMControlAlias = &fNormalizedRPMControlUpscaled;
for( int n = 0; n < iNumSamples; n++ )
{
// Low-pass input
fNormalizedRPMControlUpscaled = fNormalizedRPMControl;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMControlAlias );
fNormalizedRPMControlUpscaled = *pfNormalizedRPMControlAlias;
/* Update DSP network coefficients */
float& fCurrentSample( pfOutputBuffer[ n ] );
// Update turbine frequencies
assert( m_vfTurbineBaseModeFrequencies.size() == oCtxAudio.vfTurbineModeFrequencies.size() );
for( int i = 0; i < m_vfTurbineBaseModeFrequencies.size(); i++ )
oCtxAudio.vfTurbineModeFrequencies[ i ] = m_vfTurbineBaseModeFrequencies[ i ] * fNormalizedRPMControlUpscaled;
// Update forced flame filters
const float fCenterFrequency2 = fNormalizedRPMControlUpscaled * fNormalizedRPMControlUpscaled * 150.0f;
m_oForcedFlameBP2.setup( g_iIIRFilterOrder, float( m_dSampleRate ), fCenterFrequency2, fCenterFrequency2 / 1.0f );
const float fCenterFrequency3 = fNormalizedRPMControlUpscaled * 12000.0f;
const float fQ3 = 0.6f;
m_oForcedFlameBP3.setup( g_iIIRFilterOrder, float( m_dSampleRate ), fCenterFrequency3, fCenterFrequency3 / fQ3 );
/* Process DSP */
// Forced flame
m_fTempSample = oRNG.GenerateFloat( -1.0f, 1.0f ); // noise~
m_oForcedFlameBP1.process( 1, &pfTempSampleAlias ); // bd~ 8000 0.5
m_oForcedFlameBP2.process( 1, &pfTempSampleAlias ); // vcf~ 0 1 (real part only = band pass)
m_oForcedFlameHP1.process( 1, &pfTempSampleAlias ); // hip~ 120
m_fTempSample *= 120.0f; // *~ 120
m_fTempSample = ( m_fTempSample < -1.0f ) ? -1 : ( ( m_fTempSample > 1 ) ? 1 : m_fTempSample ); // clip~ -1 1
m_oForcedFlameBP3.process( 1, &pfTempSampleAlias ); // vcf~ 0 0.6 (real part only = band pass)
fCurrentSample = m_fTempSample; // override buffer
/* Turbine (adds to output buffer) */
m_fTempSample = 0.0f;
assert( oCtxAudio.vfTurbineModeFrequencies.size() == oCtxAudio.vfTurbineModePhaseShift.size() );
assert( oCtxAudio.vfTurbineModeFrequencies.size() == m_vfTurbineBaseModeAmplitudes.size() );
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
const float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
float fAmplitude( m_vfTurbineBaseModeAmplitudes[ i ] );
const float t = fmodf( ITAConstants::TWO_PI_F_L * fFrequency * float( n ) / float( m_dSampleRate ), ITAConstants::TWO_PI_F );
// [MANUAL MODIFICATION] sine-sqared roll-in of amplitude for the frequencies, otherwise very synthetic
if( fNormalizedRPMControlUpscaled < 0.1 )
fAmplitude *= sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f ) * sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f );
m_fTempSample += fAmplitude * sin( t + fPhaseShift ); // all osc~ and *~
}
m_fTempSample = ( m_fTempSample < -0.9f ) ? -0.9f : ( ( m_fTempSample > 0.9f ) ? 0.9f : m_fTempSample ); // clip~ -0.9 0.9
/* Jet engine */
m_fTempSample *= 0.5f; // *~ 0.5
const float fManualBalance = 0.1f / 0.5f; // [MANUAL MODIFICATION] not included in patch
m_fTempSample = fCurrentSample + fManualBalance * m_fTempSample; // combine turbine and flame from jet engine patch (with a manual balance that sound better)
m_oJetEngineLP1.process( 1, &pfTempSampleAlias ); // ~lop 11000
fCurrentSample = 0.2f * m_fTempSample; // *~0.2
}
/* Update phases for turbine */
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
fPhaseShift = fmodf( ITAConstants::TWO_PI_F_L * fFrequency / float( m_dSampleRate ) * float( iNumSamples ) + fPhaseShift, ITAConstants::TWO_PI_F );
}
}
......@@ -3,117 +3,26 @@
#include <ITASampleBuffer.h>
#include <ITAStopWatch.h>
#include <ITAStringUtils.h>
#include <ITAFiniteImpulseResponse.h>
#include <ITAMultichannelFiniteImpulseResponse.h>
#include <VistaTools/VistaRandomNumberGenerator.h>
#include <ITADSP/PD/JetEngine.h>
#include <DspFilters/Butterworth.h>
#include <atomic>
#include <algorithm>
#include <atomic>
#include <iostream>
#include "ITAFiniteImpulseResponse.h"
#include "ITAMultichannelFiniteImpulseResponse.h"
using namespace std;
using namespace Dsp;
const double g_dSampleRate = 44100;
const int g_iOutputLengthSamples = 20 * int( g_dSampleRate );
int iTotalSamples = 0;
const int g_iIIRFilterOrder = 1;
const int g_iChannels = 1;
/**
* Manual implementation of the "jet engine" pure data patch (combining "forced flame" and "turbine") from Andy Farnell's Designing Sound
*/
namespace ITADSP
{
namespace PD
{
class CJetEngine
{
public:
inline CJetEngine( float fRPMInit = 1000.0f, bool bColdStart = true )
{
m_vfTurbineBaseModeFrequencies = { 3097.0f, 4495.0f, 5588.0f, 7414.0f, 11000.0f };
m_vfTurbineBaseModeAmplitudes = { 0.25f, 0.25f, 1.0f, 0.4f, 0.4f };
m_vfRPMRange = { 1000.0f, 5000.0f };
const float fCenterFrequency1 = 8000.0f;
const float fQ1 = 0.5f;
m_oForcedFlameBP1.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCenterFrequency1, fCenterFrequency1 / fQ1 );
const float fCutoffFrequency1 = 120.0f;
m_oForcedFlameHP1.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCutoffFrequency1 );
const float fLPCutoffFrequency1 = 11000.0f;
m_oJetEngineLP1.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fLPCutoffFrequency1 );
const float fLPCutoffFrequency2 = 0.1f;
m_oJetEngineLP2.setup( 1, float( g_dSampleRate ), fLPCutoffFrequency2 );
oCtxAudio.vfTurbineModePhaseShift = { 0, 0, 0, 0, 0 };
oCtxAudio.vfTurbineModeFrequencies = { 0, 0, 0, 0, 0 };
SetRPM( fRPMInit );
if( !bColdStart )
{
// Run input low-pass for 3 secs (until 0.2 Hz LP died away on unit step function) to akkumulate
// of the default input, otherwise there is a run-up from zero to target RPM value
const float fNormalizedRPMInit = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMUpsample;
float* pfNormalizedRPMUpsampleAlias = &fNormalizedRPMUpsample;
for( int n = 0; n < ( int ) g_dSampleRate * 3; n++ )
{
fNormalizedRPMUpsample = fNormalizedRPMInit;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMUpsampleAlias );
}
}
};
inline virtual ~CJetEngine()
{
};
virtual void SetRPM( float fNewRPM );
virtual void Process( float* pfOutputBuffer, int iNumSamples );
protected:
// We'll protect data exchange between user and audio thread
// with atomics, but a request, swap & acknowledge approach with
// mutex / locking would be more elegant
struct CAudioContext
{
std::atomic< float > vfRPM;
std::vector< float > vfTurbineModeFrequencies;
std::vector< float > vfTurbineModePhaseShift;
} oCtxAudio;
float GetRPMNormalized( float fRPM ) const;
private:
VistaRandomNumberGenerator oRNG;
SimpleFilter< Butterworth::BandPass< g_iIIRFilterOrder >, g_iChannels > m_oForcedFlameBP1, m_oForcedFlameBP2, m_oForcedFlameBP3;
SimpleFilter< Butterworth::HighPass< g_iIIRFilterOrder >, g_iChannels > m_oForcedFlameHP1;
SimpleFilter< Butterworth::LowPass< 1 >, g_iChannels > m_oJetEngineLP1, m_oJetEngineLP2;
std::vector< float > m_vfTurbineBaseModeFrequencies, m_vfTurbineBaseModeAmplitudes, m_vfRPMRange;
float m_fTempSample;
};
}
}
int main( int, char** )
{
ITASampleBuffer oOutputBuffer( g_iOutputLengthSamples );
vector<float > vfRPMs = { 1000.f, 4000.0f, 2000.0f };
ITADSP::PD::CJetEngine oPatch( vfRPMs[ 0 ] );
bool bColdStart = true;
ITADSP::PD::CJetEngine oPatch( g_dSampleRate, vfRPMs[ 0 ], bColdStart );
int iTimeSeriesLeg = g_iOutputLengthSamples / ( int ) vfRPMs.size();
......@@ -133,122 +42,3 @@ int main( int, char** )
return 255;
}
void ITADSP::PD::CJetEngine::SetRPM( float fRPM )
{
oCtxAudio.vfRPM = fRPM; // Direct write-through, no lock & swap implemented here
}
float ITADSP::PD::CJetEngine::GetRPMNormalized( float fRPM ) const
{
// Validate input range
assert( m_vfRPMRange.size() == 2 );
const float fValidRPM = std::min( std::max( m_vfRPMRange[ 0 ], fRPM ), m_vfRPMRange[ 1 ] );
const float fNormalizedRPM = std::max( 0.1f, fValidRPM / m_vfRPMRange[ 1 ] );
return fNormalizedRPM;
}
void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
{
/* Validate, normalize and up-sampled RPM and apply low-pass filter */
// Processed DSP coefficient with normalized input
const float fNormalizedRPMControl = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMControlUpscaled;
// We need this ref pointer for Dsp functions, but its ugly we dont use it apart from that
float* pfTempSampleAlias = &m_fTempSample;
float* pfNormalizedRPMControlAlias = &fNormalizedRPMControlUpscaled;
for( int n = 0; n < iNumSamples; n++ )
{
// Low-pass input
fNormalizedRPMControlUpscaled = fNormalizedRPMControl;
m_oJetEngineLP2.process( 1, &pfNormalizedRPMControlAlias );
fNormalizedRPMControlUpscaled = *pfNormalizedRPMControlAlias;
/* Update DSP network coefficients */
float& fCurrentSample( pfOutputBuffer[ n ] );
// Update turbine frequencies
assert( m_vfTurbineBaseModeFrequencies.size() == oCtxAudio.vfTurbineModeFrequencies.size() );
for( int i = 0; i < m_vfTurbineBaseModeFrequencies.size(); i++ )
oCtxAudio.vfTurbineModeFrequencies[ i ] = m_vfTurbineBaseModeFrequencies[ i ] * fNormalizedRPMControlUpscaled;
// Update forced flame filters
const float fCenterFrequency2 = fNormalizedRPMControlUpscaled * fNormalizedRPMControlUpscaled * 150.0f;
m_oForcedFlameBP2.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCenterFrequency2, fCenterFrequency2 / 1.0f );
const float fCenterFrequency3 = fNormalizedRPMControlUpscaled * 12000.0f;
const float fQ3 = 0.6f;
m_oForcedFlameBP3.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCenterFrequency3, fCenterFrequency3 / fQ3 );
/* Process DSP */
// Forced flame
m_fTempSample = oRNG.GenerateFloat( -1.0f, 1.0f ); // noise~
m_oForcedFlameBP1.process( 1, &pfTempSampleAlias ); // bd~ 8000 0.5
m_oForcedFlameBP2.process( 1, &pfTempSampleAlias ); // vcf~ 0 1 (real part only = band pass)
m_oForcedFlameHP1.process( 1, &pfTempSampleAlias ); // hip~ 120
m_fTempSample *= 120.0f; // *~ 120
m_fTempSample = ( m_fTempSample < -1.0f ) ? -1 : ( ( m_fTempSample > 1 ) ? 1 : m_fTempSample ); // clip~ -1 1
m_oForcedFlameBP3.process( 1, &pfTempSampleAlias ); // vcf~ 0 0.6 (real part only = band pass)
fCurrentSample = m_fTempSample; // override buffer
/* Turbine (adds to output buffer) */
m_fTempSample = 0.0f;
assert( oCtxAudio.vfTurbineModeFrequencies.size() == oCtxAudio.vfTurbineModePhaseShift.size() );
assert( oCtxAudio.vfTurbineModeFrequencies.size() == m_vfTurbineBaseModeAmplitudes.size() );
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
const float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
float fAmplitude( m_vfTurbineBaseModeAmplitudes[ i ] );
const float t = fmodf( ITAConstants::TWO_PI_F_L * fFrequency * float( n ) / float( g_dSampleRate ), ITAConstants::TWO_PI_F );
// [MANUAL MODIFICATION] sine-sqared roll-in of amplitude for the frequencies, otherwise very synthetic
if( fNormalizedRPMControlUpscaled < 0.1 )
fAmplitude *= sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f ) * sin( ITAConstants::HALF_PI_F * fNormalizedRPMControlUpscaled / 0.1f );
m_fTempSample += fAmplitude * sin( t + fPhaseShift ); // all osc~ and *~
}
m_fTempSample = ( m_fTempSample < -0.9f ) ? -0.9f : ( ( m_fTempSample > 0.9f ) ? 0.9f : m_fTempSample ); // clip~ -0.9 0.9
/* Jet engine */
m_fTempSample *= 0.5f; // *~ 0.5
const float fManualBalance = 0.1f / 0.5f; // [MANUAL MODIFICATION] not included in patch
m_fTempSample = fCurrentSample + fManualBalance * m_fTempSample; // combine turbine and flame from jet engine patch (with a manual balance that sound better)
m_oJetEngineLP1.process( 1, &pfTempSampleAlias ); // ~lop 11000
fCurrentSample = 0.2f * m_fTempSample; // *~0.2
}
/* Update phases for turbine */
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
fPhaseShift = fmodf( ITAConstants::TWO_PI_F_L * fFrequency / float( g_dSampleRate ) * float( iNumSamples ) + fPhaseShift, ITAConstants::TWO_PI_F );
}
}
%%
jet_engine_in_params = ita_read( 'ITADSP_pd_jet_engine_input_param.wav' );
%jet_engine_in_params.pt
%%
jet_engine = ita_read( 'ITADSP_pd_jet_engine_out.wav' );
jet_engine.plot_spectrogram
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment