Skip to content

Commit

Permalink
Use log2 approximation to increase performance. Part of #2004
Browse files Browse the repository at this point in the history
  • Loading branch information
@srcejon
srcejon committed Mar 4, 2024
1 parent d48e962 commit 0b372d9
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 16 deletions.
78 changes: 62 additions & 16 deletions sdrbase/dsp/spectrumvis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
// Copyright (C) 2015-2023 Edouard Griffiths, F4EXB <[email protected]> //
// Copyright (C) 2022 Jiří Pinkava <[email protected]> //
// Copyright (C) 2023 Arne Jünemann <[email protected]> //
// Copyright (C) 2023 Vladimir Pleskonjic //
// //
// This program is free software; you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
Expand Down Expand Up @@ -135,7 +136,7 @@ void SpectrumVis::feed(const Complex *begin, unsigned int length)

v = c.real() * c.real() + c.imag() * c.imag();
m_psd[i] = v/m_powFFTDiv;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i] = v;
}

Expand Down Expand Up @@ -176,7 +177,7 @@ void SpectrumVis::feed(const Complex *begin, unsigned int length)
v = c.real() * c.real() + c.imag() * c.imag();
v = m_movingAverage.storeAndGetAvg(v, i);
m_psd[i] = v/m_powFFTDiv;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i] = v;
}

Expand Down Expand Up @@ -224,7 +225,7 @@ void SpectrumVis::feed(const Complex *begin, unsigned int length)
if (m_fixedAverage.storeAndGetAvg(avg, v, i))
{
m_psd[i] = avg/m_powFFTDiv;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2f(avg) + m_ofs;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2fapprox(avg) + m_ofs;
m_powerSpectrum[i] = avg;
}
}
Expand Down Expand Up @@ -275,7 +276,7 @@ void SpectrumVis::feed(const Complex *begin, unsigned int length)
if (m_max.storeAndGetMax(max, v, i))
{
m_psd[i] = max/m_powFFTDiv;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2f(max) + m_ofs;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2fapprox(max) + m_ofs;
m_powerSpectrum[i] = max;
}
}
Expand Down Expand Up @@ -450,7 +451,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
v = c.real() * c.real() + c.imag() * c.imag();
m_psd[i] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i * 2] = v;
m_powerSpectrum[i * 2 + 1] = v;
}
Expand All @@ -463,14 +464,14 @@ void SpectrumVis::processFFT(bool positiveOnly)
v = c.real() * c.real() + c.imag() * c.imag();
m_psd[i] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i] = v;

c = fftOut[i];
v = c.real() * c.real() + c.imag() * c.imag();
m_psd[i + halfSize] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i + halfSize] = v;
}
}
Expand Down Expand Up @@ -512,7 +513,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
v = m_movingAverage.storeAndGetAvg(v, i);
m_psd[i] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i * 2] = v;
m_powerSpectrum[i * 2 + 1] = v;
}
Expand All @@ -526,15 +527,15 @@ void SpectrumVis::processFFT(bool positiveOnly)
v = m_movingAverage.storeAndGetAvg(v, i+halfSize);
m_psd[i] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i] = v;

c = fftOut[i];
v = c.real() * c.real() + c.imag() * c.imag();
v = m_movingAverage.storeAndGetAvg(v, i);
m_psd[i + halfSize] = v/m_powFFTDiv;
m_specMax = v > m_specMax ? v : m_specMax;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2f(v) + m_ofs;
v = m_settings.m_linear ? v/m_powFFTDiv : m_mult * log2fapprox(v) + m_ofs;
m_powerSpectrum[i + halfSize] = v;
}
}
Expand Down Expand Up @@ -582,7 +583,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i] = avg/m_powFFTDiv;
specMax = avg > specMax ? avg : specMax;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2f(avg) + m_ofs;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2fapprox(avg) + m_ofs;
m_powerSpectrum[i * 2] = avg;
m_powerSpectrum[i * 2 + 1] = avg;
}
Expand All @@ -600,7 +601,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i] = avg/m_powFFTDiv;
specMax = avg > specMax ? avg : specMax;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2f(avg) + m_ofs;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2fapprox(avg) + m_ofs;
m_powerSpectrum[i] = avg;
}

Expand All @@ -612,7 +613,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i + halfSize] = avg/m_powFFTDiv;
specMax = avg > specMax ? avg : specMax;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2f(avg) + m_ofs;
avg = m_settings.m_linear ? avg/m_powFFTDiv : m_mult * log2fapprox(avg) + m_ofs;
m_powerSpectrum[i + halfSize] = avg;
}
}
Expand Down Expand Up @@ -665,7 +666,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i] = max/m_powFFTDiv;
specMax = max > specMax ? max : specMax;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2f(max) + m_ofs;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2fapprox(max) + m_ofs;
m_powerSpectrum[i * 2] = max;
m_powerSpectrum[i * 2 + 1] = max;
}
Expand All @@ -683,7 +684,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i] = max/m_powFFTDiv;
specMax = max > specMax ? max : specMax;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2f(max) + m_ofs;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2fapprox(max) + m_ofs;
m_powerSpectrum[i] = max;
}

Expand All @@ -695,7 +696,7 @@ void SpectrumVis::processFFT(bool positiveOnly)
{
m_psd[i + halfSize] = max/m_powFFTDiv;
specMax = max > specMax ? max : specMax;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2f(max) + m_ofs;
max = m_settings.m_linear ? max/m_powFFTDiv : m_mult * log2fapprox(max) + m_ofs;
m_powerSpectrum[i + halfSize] = max;
}
}
Expand Down Expand Up @@ -1084,3 +1085,48 @@ void SpectrumVis::webapiUpdateSpectrumSettings(

settings.updateFrom(prefixedKeys, &response);
}

// To calculate power, the usual equation:
// 10*log10(V1/V2), where V2=fftSize^2
// is calculated using log2 instead, with:
// ofs=20.0f * log10f(1.0f / fftSize)
// mult=(10.0f / log2(10.0f))
// dB = m_mult * log2f(v) + m_ofs
// However, while the gcc version of log2f is twice as fast as log10f,
// MSVC version is 6x slower.
// Also, we don't need full accuracy of log2f for calculating the power for the spectrum,
// so we can use the following approximation to get a good speed-up for both compilers:
// https://www.vplesko.com/posts/replacing_log2f.html
// https://www.vplesko.com/assets/replacing_log2f/main.c.txt
float SpectrumVis::log2fapprox(float x)
{
// IEEE 754 representation constants.
const int32_t mantissaLen = 23;
const int32_t mantissaMask = (1 << mantissaLen) - 1;
const int32_t baseExponent = -127;

// Reinterpret x as int in a standard compliant way.
int32_t xi;
memcpy(&xi, &x, sizeof(xi));

// Calculate exponent of x.
float e = (float)((xi >> mantissaLen) + baseExponent);

// Calculate mantissa of x. It will be in range [1, 2).
float m;
int32_t mxi = (xi & mantissaMask) | ((-baseExponent) << mantissaLen);
memcpy(&m, &mxi, sizeof(m));

// Use Remez algorithm-generated approximation polynomial
// for log2(a) where a is in range [1, 2].
float l = 0.15824871f;
l = l * m + -1.051875f;
l = l * m + 3.0478842f;
l = l * m + -2.1536207f;

// Add exponent to the calculation.
// Final log is log2(m*2^e)=log2(m)+e.
l += e;

return l;
}
1 change: 1 addition & 0 deletions sdrbase/dsp/spectrumvis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,7 @@ class SDRBASE_API SpectrumVis : public QObject, public BasebandSampleSink {
void handleScalef(Real scalef);
void handleWSOpenClose(bool openClose);
void handleConfigureWSSpectrum(const QString& address, uint16_t port);
float log2fapprox(float x);

static void webapiFormatSpectrumSettings(SWGSDRangel::SWGGLSpectrum& response, const SpectrumSettings& settings);
static void webapiUpdateSpectrumSettings(
Expand Down

0 comments on commit 0b372d9

Please sign in to comment.