ITAThirdOctaveFIRFilterGenerator.cpp 4.72 KB
Newer Older
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
1
2
3
4
5
6
7
8
9
10
#include <ITAThirdOctaveFIRFilterGenerator.h>

#include <ITAConstants.h>
#include <ITAFastMath.h>
#include <ITANumericUtils.h>
#include <ITAStringUtils.h>
#include <ITAThirdOctaveMagnitudeSpectrum.h>

#include <spline.h>

11
12
using namespace ITABase;

Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
13
14
15
16
17
18
19
20
21
22
CITAThirdOctaveFIRFilterGenerator::CITAThirdOctaveFIRFilterGenerator( const double dSampleRate, const int iFilterLength )
	: m_dSamplerate( dSampleRate )
	, m_iFilterLength( iFilterLength )
	, m_ypp( nullptr )
	, m_pfInputFreqs( nullptr )
	, m_pfInputData( nullptr )
	, m_pfBuf1( nullptr )
	, m_pfBuf2( nullptr )
	, m_bWindow( false )
{
23
	m_iInputFreqs = CThirdOctaveGainMagnitudeSpectrum::GetNumBands() + 2;
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
24
25
	m_pfInputFreqs = fm_falloc( m_iInputFreqs, true );
	m_pfInputFreqs[ 0 ] = 0;	// Left margin
26
27
	for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
		m_pfInputFreqs[ i + 1 ] = CThirdOctaveGainMagnitudeSpectrum::GetCenterFrequencies()[ i ];
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
	m_pfInputFreqs[ m_iInputFreqs - 1 ] = ( float ) dSampleRate / 2;	// Right margin: Nyquist frequency

	m_pfInputData = fm_falloc( m_iInputFreqs, true );

	// DFT frequency bandwidth
	m_fDeltaF = ( float ) dSampleRate / ( float ) iFilterLength;

	// Number of symetric DFT coefficients;
	m_iDFTCoeffs = iFilterLength / 2 + 1;

	m_pfBuf1 = fm_falloc( 2 * m_iDFTCoeffs, false );
	m_pfBuf2 = fm_falloc( iFilterLength, false );
	m_pfWindow = fm_falloc( iFilterLength, false );

	// Windowing function (Hann window)
	float c = 2 * ITAConstants::PI_F / ( float ) ( m_iFilterLength - 1 );
	for( int i = 0; i < m_iFilterLength; i++ )
	{
46
		m_pfWindow[ i ] = 0.5F * ( 1 - cos( c * i ) );
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
	}

	m_ifft.plan( ITAFFT::IFFT_C2R, iFilterLength, m_pfBuf1, m_pfBuf2 );

	m_sDumpFilename = "interpolated_magnitudes.csv";
}

CITAThirdOctaveFIRFilterGenerator::~CITAThirdOctaveFIRFilterGenerator()
{
	fm_free( m_pfInputFreqs );
	fm_free( m_pfInputData );
	fm_free( m_pfBuf1 );
	fm_free( m_pfBuf2 );
	fm_free( m_pfWindow );
}

int CITAThirdOctaveFIRFilterGenerator::GetFilterLength() const
{
	return m_iFilterLength;
}
67
68
69
70
71
72

double CITAThirdOctaveFIRFilterGenerator::GetSampleRate() const
{
	return m_dSamplerate;
}

Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
int CITAThirdOctaveFIRFilterGenerator::GetLatency() const
{
	// Latency = Half DFT period (ceil)
	return uprdiv( m_iFilterLength, 2 );
}

double CITAThirdOctaveFIRFilterGenerator::GetAverageRuntime() const
{
	return m_sw.mean();
}

void CITAThirdOctaveFIRFilterGenerator::SetDumpFilename( const std::string& sFilename )
{
	m_sDumpFilename = sFilename;
}

89
void CITAThirdOctaveFIRFilterGenerator::GenerateFilter( const ITABase::CThirdOctaveGainMagnitudeSpectrum& oTOGainMagnitudes, float* pfFilterCoeffs, bool bMinimumPhase /*=false*/ )
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
90
91
92
{
	m_sw.start();

93
	if( oTOGainMagnitudes.IsZero() )
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
94
	{
95
96
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = 0.0f;
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
97
98
99
		return;
	}

100
	if( oTOGainMagnitudes.IsIdentity() )
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
101
	{
102
103
104
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = 0.0f;
		pfFilterCoeffs[ int( m_iFilterLength / 2 ) ] = 1.0f;
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
105
106
		return;
	}
107

Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
108
109
	// 1st step: Interpolate the magnitudes

110
111
	for( int i = 0; i < CThirdOctaveGainMagnitudeSpectrum::GetNumBands(); i++ )
		m_pfInputData[ 1 + i ] = float( oTOGainMagnitudes[ i ] );
112
113
114
115

	// Bounndaries must be defined in a meaningful way.
	m_pfInputData[ 0 ] = m_pfInputData[ 1 ];
	m_pfInputData[ m_iInputFreqs - 1 ] = 0.0f;
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
116
117

	// Initialize cubic spline interpolation
118
	m_ypp = spline_cubic_set( m_iInputFreqs,
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
119
120
121
122
123
		m_pfInputFreqs,
		m_pfInputData,
		1, // Left boundary condition => 1st derivative m=0
		0,
		1, // Right boundary condition => 1st derivative m=0
124
		0 );
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
125
	float fDummy;
126
	const float fScale = 1 / ( float ) m_iFilterLength;
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
127
128

	// No DC offset, ever!
129
130
	m_pfBuf1[ 0 ] = 0;
	m_pfBuf1[ 1 ] = 0;
131
132


133
134
	if( bMinimumPhase ) {
		for( int i = 1; i < m_iDFTCoeffs; i++ )
135
		{
136
			float x = spline_cubic_val( m_iInputFreqs,
137
138
139
140
141
				m_pfInputFreqs,
				i*m_fDeltaF,
				m_pfInputData,
				m_ypp,
				&fDummy,
142
				&fDummy );
143
144

			// Phase-shift by half the FFT-period
145
146
			m_pfBuf1[ 2 * i ] = pow( x * fScale, 2 ) * m_iFilterLength; //minimum phase
			m_pfBuf1[ 2 * i + 1 ] = 0;
147
148
149
		}
	}
	else {
150
		for( int i = 1; i < m_iDFTCoeffs; i++ )
151
		{
152
			float x = spline_cubic_val( m_iInputFreqs,
153
154
155
156
157
				m_pfInputFreqs,
				i*m_fDeltaF,
				m_pfInputData,
				m_ypp,
				&fDummy,
158
				&fDummy );
159
160

			// Phase-shift by half the FFT-period: Negate all odd DFT coefficients
161
162
			m_pfBuf1[ 2 * i ] = ( ( i % 2 ) == 0 ) ? x * fScale : -x * fScale;
			m_pfBuf1[ 2 * i + 1 ] = 0;
163
		}
Dipl.-Ing. Jonas Stienen's avatar
Dipl.-Ing. Jonas Stienen committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
	}

	// 2nd step: Convert into time-domain (out-of-place C2R-IFFT)
	m_ifft.execute();

	// 3rd (optional) step: Hann window in the time-domain (optional)
	if( m_bWindow )
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = m_pfBuf2[ i ] * m_pfWindow[ i ];
	}
	else
	{
		for( int i = 0; i < m_iFilterLength; i++ )
			pfFilterCoeffs[ i ] = m_pfBuf2[ i ];
	}

	// @todo: Minimum-phase?

	m_sw.stop();
}