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

Fixing problems with jet engine noise and figuring out problems with single...

Fixing problems with jet engine noise and figuring out problems with single zero pole low pass for RPM input
parent a04aa5a6
......@@ -48,6 +48,7 @@ void tri_band();
void downcast_test();
void octave_bands();
void live_coeff_change_test();
void very_low_low_pass();
int main( int, char** )
{
......@@ -63,6 +64,8 @@ int main( int, char** )
live_coeff_change_test();
very_low_low_pass();
return 255;
}
......@@ -235,3 +238,46 @@ void live_coeff_change_test()
sbCoeffChangeTest.Store( "NativeDspFiltersTest_LiveCoeffChangeTest.wav" );
}
void very_low_low_pass()
{
float fVeryLowLowPassFrequencyCuttoff = 0.2f; // 0.2 Hz
SimpleFilter< Butterworth::LowPass< 1 >, 1 > oVLLPF;
oVLLPF.setup( 1, fSampleRate, fVeryLowLowPassFrequencyCuttoff );
CFiniteImpulseResponse oSF( iFilterLength * 30, fSampleRate );
oSF.SetUnitStepFunction( 0.1f );
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 2 * iFilterLength + n ] = 0.2f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 4 * iFilterLength + n ] = 0.3f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 6 * iFilterLength + n ] = 0.4f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 8 * iFilterLength + n ] = 0.5f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 10 * iFilterLength + n ] = 1.0f;
for( int n = 0; n < 4 * iFilterLength; n++ )
oSF[ 2 * 12 * iFilterLength + n ] = 0.0f;
oSF.Store( "NativeDspFiltersTest_SF_VeryLowLowPass_in.wav" );
for( size_t n = 0; n < oSF.GetLength(); n++ )
{
auto pf = oSF.GetData() + n;
oVLLPF.process( 1, &pf );
}
oSF.Store( "NativeDspFiltersTest_SF_VeryLowLowPass_out.wav" );
oVLLPF.reset();
CFiniteImpulseResponse oIR( iFilterLength * 30, fSampleRate );
oIR.SetDirac();
for( size_t n = 0; n < oIR.GetLength(); n++ )
{
auto pf = oIR.GetData() + n;
oVLLPF.process( 1, &pf);
}
oIR.Store( "NativeDspFiltersTest_IR_VeryLowLowPass.wav" );
}
......@@ -49,3 +49,16 @@ IR_Octaves.pf
LiveCoeffChangeTest = ita_read( 'NativeDspFiltersTest_LiveCoeffChangeTest.wav' );
LiveCoeffChangeTest.pt
LiveCoeffChangeTest.ch(1).play
%% Very low low pass
vllp_ir = ita_read( 'NativeDspFiltersTest_IR_VeryLowLowPass.wav' );
vllp_sf_in = ita_read( 'NativeDspFiltersTest_SF_VeryLowLowPass_in.wav' );
vllp_sf_out = ita_read( 'NativeDspFiltersTest_SF_VeryLowLowPass_out.wav' );
vllp = ita_merge( vllp_ir, vllp_sf_in, vllp_sf_out );
%vllp.channelUnits = 'energy';
vllp.channelNames = { 'Impulse response', 'Step function [in]', 'Step function [out]' };
title 'Very low pass filter (2 Hz) first order zero pole'
vllp.pt
......@@ -8,15 +8,19 @@
#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 = 10 * int( g_dSampleRate );
const int g_iOutputLengthSamples = 20 * int( g_dSampleRate );
int iTotalSamples = 0;
const int g_iIIRFilterOrder = 1;
const int g_iChannels = 1;
......@@ -31,13 +35,12 @@ namespace ITADSP
class CJetEngine
{
public:
inline CJetEngine( float RPMInit = 500.0f )
inline CJetEngine( float fRPMInit = 1000.0f, bool bColdStart = true )
{
m_vfTurbineModeFrequencies = { 3097.0f, 4495.0f, 5588.0f, 7414.0f, 11000.0f };
m_vfTurbineModeAmplitudes = { 0.25f, 0.25f, 1.0f, 0.4f, 0.4f };
m_vfTurbineModePhaseShift = { 0, 0, 0, 0, 0 };
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 = { 500.0f, 6000.0f };
m_vfRPMRange = { 1000.0f, 5000.0f };
const float fCenterFrequency1 = 8000.0f;
const float fQ1 = 0.5f;
......@@ -49,28 +52,64 @@ namespace ITADSP
const float fLPCutoffFrequency1 = 11000.0f;
m_oJetEngineLP1.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fLPCutoffFrequency1 );
SetRPM( RPMInit );
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 );
}
}
oInputParameter.Init( 4, g_iOutputLengthSamples, g_dSampleRate, true );
};
inline virtual ~CJetEngine() {};
inline virtual ~CJetEngine()
{
oInputParameter.StoreToFile( "ITADSP_pd_jet_engine_input_param.wav" );
};
virtual void SetRPM( float fNewRPM );
virtual void Process( float* pfOutputBuffer, int iNumSamples );
protected:
void UpdateForcedFlame( float fRPM );
void UpdateTurbine( float fRPM );
// 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< g_iIIRFilterOrder >, g_iChannels > m_oJetEngineLP1;
SimpleFilter< Butterworth::LowPass< 1 >, g_iChannels > m_oJetEngineLP1, m_oJetEngineLP2;
std::vector< float > m_vfTurbineModeFrequencies, m_vfRPMRange, m_vfTurbineModePhaseShift, m_vfTurbineModeAmplitudes;
std::vector< float > m_vfTurbineBaseModeFrequencies, m_vfTurbineBaseModeAmplitudes, m_vfRPMRange;
float m_fTempSample;
ITABase::CMultichannelFiniteImpulseResponse oInputParameter;
};
}
}
......@@ -79,16 +118,16 @@ int main( int, char** )
{
ITASampleBuffer oOutputBuffer( g_iOutputLengthSamples );
vector<float > vfRPMs = { 700.f, 800.0f, 600.0f, 1100.0f, 2500.0f, 4000.0f, 4200.0f, 4000.0f, 3000.0f, 2800.0f };
vector<float > vfRPMs = { 1000.f, 4000.0f, 2000.0f };
ITADSP::PD::CJetEngine oPatch( vfRPMs[ 0 ] );
int iTimeSeriesLeg = g_iOutputLengthSamples / (int) vfRPMs.size();
int iTotalSamples = 0;
int iTimeSeriesLeg = g_iOutputLengthSamples / ( int ) vfRPMs.size();
for( auto m = 0; m < vfRPMs.size(); m++ )
{
oPatch.SetRPM( vfRPMs[ m ] );
int iProcessSamples = std::min( iTimeSeriesLeg, oOutputBuffer.GetLength() - m * iTimeSeriesLeg );
int iProcessSamples = std::min( iTimeSeriesLeg, g_iOutputLengthSamples - m * iTimeSeriesLeg );
oPatch.Process( oOutputBuffer.GetData() + iTotalSamples, iProcessSamples );
iTotalSamples += iProcessSamples;
}
......@@ -103,46 +142,57 @@ int main( int, char** )
void ITADSP::PD::CJetEngine::SetRPM( float fRPM )
{
UpdateForcedFlame( fRPM );
UpdateTurbine( fRPM );
oCtxAudio.vfRPM = fRPM; // Direct write-through, no lock & swap implemented here
}
void ITADSP::PD::CJetEngine::UpdateForcedFlame( float fRPM )
float ITADSP::PD::CJetEngine::GetRPMNormalized( float fRPM ) const
{
// Validate input range
assert( m_vfRPMRange.size() == 2 );
float fValidRPM = std::min( std::max( m_vfRPMRange[ 0 ], fRPM ), m_vfRPMRange[ 1 ] );
// Normalize input
const float fNormalizedFlameControl = std::max( 0.1f, fValidRPM / m_vfRPMRange[ 1 ] );
const float fCenterFrequency2 = fNormalizedFlameControl * fNormalizedFlameControl * 150.0f;
m_oForcedFlameBP2.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCenterFrequency2, fCenterFrequency2 / 1.0f );
const float fCenterFrequency3 = fNormalizedFlameControl * 12000.0f;
const float fQ3 = 0.6f;
m_oForcedFlameBP3.setup( g_iIIRFilterOrder, float( g_dSampleRate ), fCenterFrequency3, fCenterFrequency3 / fQ3 );
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::UpdateTurbine( float fRPM )
void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
{
assert( m_vfRPMRange.size() == 2 );
float fValidRPM = std::min( std::max( m_vfRPMRange[ 0 ], fRPM ), m_vfRPMRange[ 1 ] );
/* Validate, normalize and up-sampled RPM and apply low-pass filter */
// Normalize input
const float fNormalizedTurbineControl = std::max( 0.1f, fValidRPM / m_vfRPMRange[ 1 ] );
// Processed DSP coefficient with normalized input
const float fNormalizedRPMControl = GetRPMNormalized( oCtxAudio.vfRPM );
float fNormalizedRPMControlUpscaled;
for( auto& v : m_vfTurbineModeFrequencies )
v *= fNormalizedTurbineControl;
}
void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
{
float* pfTempSampleAlias = &m_fTempSample; // We need this ref pointer for Dsp functions, but its ugly we dont use it apart from that
// 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~
......@@ -162,45 +212,57 @@ void ITADSP::PD::CJetEngine::Process( float* pfOutputBuffer, int iNumSamples )
fCurrentSample = m_fTempSample; // override buffer
// Turbine (adds to output buffer)
/* Turbine (adds to output buffer) */
m_fTempSample = 0.0f;
assert( m_vfTurbineModeFrequencies.size() == m_vfTurbineModePhaseShift.size() );
for( int i = 0; i < m_vfTurbineModeFrequencies.size(); i++ )
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( m_vfTurbineModeFrequencies[ i ] );
const float& fPhaseShift( m_vfTurbineModePhaseShift[ i ] );
const float& fAmplitude( m_vfTurbineModeAmplitudes[ 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 );
const int iPeriodLengthSamples = ( int ) round( g_dSampleRate / fFrequency );
const float t = fmodf( ITAConstants::TWO_PI_F_L / float( iPeriodLengthSamples ) * float( n ), ITAConstants::TWO_PI_F );
// [MANUAL WORK] sin-sqared roll-in of amplitude for the frequencies
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 *~
if( i == 0 )
{
oInputParameter[ 0 ][ iTotalSamples + n ] = m_fTempSample;
oInputParameter[ 1 ][ iTotalSamples + n ] = fFrequency;
oInputParameter[ 2 ][ iTotalSamples + n ] = fPhaseShift;
oInputParameter[ 3 ][ iTotalSamples + n ] = fNormalizedRPMControl;
}
}
m_fTempSample = ( m_fTempSample < -0.9f ) ? -0.9f : ( ( m_fTempSample > 0.9f ) ? 0.9f : m_fTempSample ); // clip~ -0.9 0.9
// Jet engine
/* Jet engine */
m_fTempSample *= 0.5f; // *~ 0.5
const float fManualBalance = 0.2f / 0.5f; // Manual modifier [NOT INCLUDED IN PD PATCH]
fCurrentSample += fManualBalance * m_fTempSample; // combine turbine and flame from jet engine patch (with a manual balance that sound better)
const float fManualBalance = 0.1f / 0.5f; // Manual modifier [NOT INCLUDED IN PD 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; // *~0.2
fCurrentSample = 0.2f * m_fTempSample; // *~0.2
}
// Update phases for turbine
for( int i = 0; i < m_vfTurbineModeFrequencies.size(); i++ )
/* Update phases for turbine */
for( int i = 0; i < oCtxAudio.vfTurbineModeFrequencies.size(); i++ )
{
const float& fFrequency( m_vfTurbineModeFrequencies[ i ] );
float& fPhaseShift( m_vfTurbineModePhaseShift[ i ] );
const float& fFrequency( oCtxAudio.vfTurbineModeFrequencies[ i ] );
float& fPhaseShift( oCtxAudio.vfTurbineModePhaseShift[ i ] );
const int iPeriodLengthSamples = ( int ) round( g_dSampleRate / fFrequency );
fPhaseShift = fmodf( ITAConstants::TWO_PI_F_L / float( iPeriodLengthSamples ) * iNumSamples + fPhaseShift, ITAConstants::TWO_PI_F );
fPhaseShift = fmodf( ITAConstants::TWO_PI_F_L * fFrequency / float( g_dSampleRate ) * float( iNumSamples ) + fPhaseShift, ITAConstants::TWO_PI_F );
}
}
jet_engine = ita_read( 'ITADSP_pd_jet_engine_out.wav' )
%%
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
jet_engine.play
\ No newline at end of file
%%
jet_engine.play
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