improved RMS calc

This commit is contained in:
Christian R. Helmrich 2020-03-22 01:00:20 +01:00
parent 27c0818160
commit 7185ac995f
4 changed files with 101 additions and 27 deletions

View File

@ -129,6 +129,7 @@ static inline uint8_t brModeAndFsToMaxSfbShort(const unsigned bitRateMode, const
return (samplingRate > 51200 ? 11 : 13) - 2 + (bitRateMode >> 2);
}
#if !SA_IMPROVED_REAL_ABS
static inline uint32_t getComplexRmsValue (const uint32_t rmsValue, const unsigned sfbGroup, const unsigned sfbIndex,
const uint8_t numSwb, const TnsData& tnsData)
{
@ -136,6 +137,7 @@ static inline uint32_t getComplexRmsValue (const uint32_t rmsValue, const unsign
return ((tnsData.numFilters > 0) && (sfbGroup == tnsData.filteredWindow) && (rmsValue <= UINT_MAX / 3) &&
(tnsData.filterLength[0] + sfbIndex >= numSwb) ? (rmsValue * 3u) >> 1 : rmsValue);
}
#endif
// ISO/IEC 23003-3, Table 75
static inline unsigned toFrameLength (const USAC_CCFL coreCoderFrameLength)
@ -374,7 +376,8 @@ unsigned ExhaleEncoder::applyTnsToWinGroup (TnsData& tnsData, SfbGroupData& grpD
memcpy (&mdctSignal[offs - MAX_PREDICTION_ORDER], m_tempIntBuf, MAX_PREDICTION_ORDER * sizeof (int32_t));
// recalculate SFB RMS in TNS range
errorValue |= m_specAnalyzer.getMeanAbsValues (mdctSignal, nullptr, nSamplesInFrame, 0, &grpSO[tnsStartSfb], __max (0, tnsMaxBands - (int) tnsStartSfb),
errorValue |= m_specAnalyzer.getMeanAbsValues (mdctSignal, nullptr /*MDST wasn't filtered*/, grpSO[grpData.sfbsPerGroup],
0 /*ci*/, &grpSO[tnsStartSfb], __max (0, tnsMaxBands - (int) tnsStartSfb),
&grpData.sfbRmsValues[tnsStartSfb + m_numSwbShort * tnsData.filteredWindow]);
}
else tnsData.filterOrder[0] = tnsData.numFilters = 0; // disable zero-length TNS filters
@ -383,10 +386,17 @@ unsigned ExhaleEncoder::applyTnsToWinGroup (TnsData& tnsData, SfbGroupData& grpD
return errorValue;
}
unsigned ExhaleEncoder::eightShortGrouping (SfbGroupData& grpData, uint16_t* const grpOffsets, int32_t* const mdctSignal)
unsigned ExhaleEncoder::eightShortGrouping (SfbGroupData& grpData, uint16_t* const grpOffsets, int32_t* const mdctSignal
#if SA_IMPROVED_REAL_ABS
, int32_t* const mdstSignal /*= nullptr*/
#endif
)
{
const unsigned nSamplesInFrame = toFrameLength (m_frameLength);
const unsigned nSamplesInShort = nSamplesInFrame >> 3;
#if SA_IMPROVED_REAL_ABS
int32_t* const tempIntBuf/*2*/ = m_timeSignals[1]; // NOTE: requires at least stereo input
#endif
unsigned grpStartLine = nSamplesInFrame;
if ((grpOffsets == nullptr) || (mdctSignal == nullptr))
@ -399,7 +409,9 @@ unsigned ExhaleEncoder::eightShortGrouping (SfbGroupData& grpData, uint16_t* con
const unsigned grpLength = grpData.windowGroupLength[gr];
uint16_t* const grpOffset = &grpOffsets[m_numSwbShort * gr];
int32_t* const grpMdctSig = &mdctSignal[grpStartLine -= nSamplesInShort * grpLength];
#if SA_IMPROVED_REAL_ABS
int32_t* const grpMdstSig = (mdstSignal != nullptr ? &mdstSignal[grpStartLine] : nullptr);
#endif
for (uint16_t b = 0; b < m_numSwbShort; b++)
{
const unsigned swbOffset = grpOffsets[b];
@ -411,13 +423,24 @@ unsigned ExhaleEncoder::eightShortGrouping (SfbGroupData& grpData, uint16_t* con
for (uint16_t w = 0; w < grpLength; w++)
{
memcpy (&m_tempIntBuf[grpOffset[b] + w * numCoeffs], &grpMdctSig[swbOffset + w * nSamplesInShort], numCoeffs * sizeof (int32_t));
#if SA_IMPROVED_REAL_ABS
if (grpMdstSig != nullptr)
{
memcpy (&tempIntBuf[grpOffset[b] + w * numCoeffs], &grpMdstSig[swbOffset + w * nSamplesInShort], numCoeffs * sizeof (int32_t));
}
#endif
}
}
grpOffset[m_numSwbShort] = uint16_t (grpStartLine + nSamplesInShort * grpLength);
} // for gr
memcpy (mdctSignal, m_tempIntBuf, nSamplesInFrame * sizeof (int32_t));
#if SA_IMPROVED_REAL_ABS
if (mdstSignal != nullptr)
{
memcpy (mdstSignal, tempIntBuf, nSamplesInFrame * sizeof (int32_t));
}
#endif
return 0; // no error
}
@ -608,10 +631,14 @@ unsigned ExhaleEncoder::psychBitAllocation () // perceptual bit-allocation via s
s = (eightShorts ? (nSamplesInFrame * grpData.windowGroupLength[gr]) >> 1 : nSamplesInFrame << 2);
for (b = 0; b < grpData.sfbsPerGroup; b++)
{
#if SA_IMPROVED_REAL_ABS
const uint32_t rmsComp = grpRms[b];
const uint32_t rmsRef9 = (coreConfig.commonWindow ? refRms[b] >> 9 : rmsComp);
#else
const uint32_t rmsComp = getComplexRmsValue (grpRms[b], gr, b, numSwbFrame, coreConfig.tnsData[ch]);
const uint32_t rmsRef9 = (!coreConfig.commonWindow ? rmsComp :
getComplexRmsValue (refRms[b], gr, b, numSwbFrame, coreConfig.tnsData[1 - ch]) >> 9);
#endif
if (rmsComp < grpRmsMin) grpRmsMin = rmsComp;
if (rmsComp >= rmsRef9 && (rmsComp < (grpStepSizes[b] >> 1))) // zero-quantized
{
@ -620,10 +647,14 @@ unsigned ExhaleEncoder::psychBitAllocation () // perceptual bit-allocation via s
}
if ((samplingRate >= 27713) && (b < maxSfbLong) && !eightShorts) // uncoded coefs
{
#if SA_IMPROVED_REAL_ABS
const uint32_t rmsComp = grpRms[b];
const uint32_t rmsRef9 = (coreConfig.commonWindow ? refRms[b] >> 9 : rmsComp);
#else
const uint32_t rmsComp = getComplexRmsValue (grpRms[b], gr, b, numSwbFrame, coreConfig.tnsData[ch]);
const uint32_t rmsRef9 = (!coreConfig.commonWindow ? rmsComp :
getComplexRmsValue (refRms[b], gr, b, numSwbFrame, coreConfig.tnsData[1 - ch]) >> 9);
#endif
if (rmsComp >= rmsRef9) // check only first SFB above max_sfb for simplification
{
s -= ((grpOff[maxSfbLong] - grpOff[b]) * reductionFactor * __min (2 * SA_EPS, rmsComp) + SA_EPS) >> 11; // / (2 * SA_EPS)
@ -978,6 +1009,8 @@ unsigned ExhaleEncoder::spectralProcessing () // complete ics_info(), calc TNS
}
else // SCE or CPE: bandwidth-to-max_sfb mapping, short-window grouping for each channel
{
coreConfig.stereoConfig = coreConfig.stereoMode = 0;
for (unsigned ch = 0; ch < nrChannels; ch++) // channel loop
{
SfbGroupData& grpData = coreConfig.groupingData[ch];
@ -999,6 +1032,8 @@ unsigned ExhaleEncoder::spectralProcessing () // complete ics_info(), calc TNS
if (samplingRate > 32000) // set max_sfb based on VBR mode and bandwidth detection
{
icsCurr.maxSfb = __min (icsCurr.maxSfb, brModeAndFsToMaxSfbLong (m_bitRateMode, samplingRate));
if (grpData.sfbsPerGroup > 49) grpData.sfbsPerGroup = 49; // for 44.1 and 48 kHz
}
while (grpSO[icsCurr.maxSfb] > __max (m_bandwidCurr[ci], m_bandwidPrev[ci])) icsCurr.maxSfb--; // BW detector
}
@ -1033,8 +1068,11 @@ unsigned ExhaleEncoder::spectralProcessing () // complete ics_info(), calc TNS
memcpy (grpData.windowGroupLength, windowGroupingTable[icsCurr.windowGrouping], NUM_WINDOW_GROUPS * sizeof (uint8_t));
#endif
while (grpSO[icsCurr.maxSfb] > __max (m_bandwidCurr[ci], m_bandwidPrev[ci])) icsCurr.maxSfb--; // not a bug!!
#if SA_IMPROVED_REAL_ABS
errorValue |= eightShortGrouping (grpData, grpSO, m_mdctSignals[ci], nChannels < 2 ? nullptr : m_mdstSignals[ci]);
#else
errorValue |= eightShortGrouping (grpData, grpSO, m_mdctSignals[ci]);
#endif
} // if EIGHT_SHORT
// compute and quantize optimal TNS coefficients, then find optimal TNS filter order
@ -1098,7 +1136,6 @@ unsigned ExhaleEncoder::spectralProcessing () // complete ics_info(), calc TNS
maxSfb0 = maxSfb1 = maxSfbSte;
}
coreConfig.commonMaxSfb = (maxSfb0 == maxSfb1); // synch
coreConfig.stereoConfig = coreConfig.stereoMode = 0;
} // if coreConfig.commonWindow
}
@ -1129,9 +1166,19 @@ unsigned ExhaleEncoder::spectralProcessing () // complete ics_info(), calc TNS
if (icsCurr.maxSfb > 0)
{
// use MCLTs for LONG but only MDCTs for SHORT windows (since MDSTs are not grouped)
// use MCLTs for LONG but only MDCTs for SHORT windows when the MDSTs aren't grouped
#if SA_IMPROVED_REAL_ABS
for (uint16_t gr = 0; gr < grpData.numWindowGroups; gr++)
{
s = m_numSwbShort * gr;
errorValue |= m_specAnalyzer.getMeanAbsValues (m_mdctSignals[ci], eightShorts && nChannels < 2 ? nullptr : m_mdstSignals[ci],
grpSO[grpData.sfbsPerGroup + s], (eightShorts ? USAC_MAX_NUM_CHANNELS : ci),
&grpSO[s], grpData.sfbsPerGroup, &grpData.sfbRmsValues[s]);
}
#else
errorValue |= m_specAnalyzer.getMeanAbsValues (m_mdctSignals[ci], eightShorts ? nullptr : m_mdstSignals[ci], nSamplesInFrame,
ci, grpSO, grpData.sfbsPerGroup * grpData.numWindowGroups, grpData.sfbRmsValues);
#endif
errorValue |= applyTnsToWinGroup (coreConfig.tnsData[ch], grpData, eightShorts, icsCurr.maxSfb, ci);
coreConfig.tnsActive |= (coreConfig.tnsData[ch].numFilters > 0); // tns_data_present
}

View File

@ -113,7 +113,11 @@ private:
// helper functions
unsigned applyTnsToWinGroup (TnsData& tnsData, SfbGroupData& grpData, const bool eightShorts, const uint8_t maxSfb,
const unsigned channelIndex);
unsigned eightShortGrouping (SfbGroupData& grpData, uint16_t* const grpOffsets, int32_t* const mdctSignal);
unsigned eightShortGrouping (SfbGroupData& grpData, uint16_t* const grpOffsets, int32_t* const mdctSignal
#if SA_IMPROVED_REAL_ABS
, int32_t* const mdstSignal = nullptr
#endif
);
unsigned getOptParCorCoeffs (const int32_t* const mdctSignal, const SfbGroupData& grpData, const uint8_t maxSfb,
const unsigned channelIndex, TnsData& tnsData, const uint8_t firstGroupIndexToTest = 0);
unsigned psychBitAllocation ();

View File

@ -1,5 +1,5 @@
/* specAnalysis.cpp - source file for class providing spectral analysis of MCLT signals
* written by C. R. Helmrich, last modified in 2019 - see License.htm for legal notices
* written by C. R. Helmrich, last modified in 2020 - see License.htm for legal notices
*
* The copyright in this software is being made available under a Modified BSD-Style License
* and comes with ABSOLUTELY NO WARRANTY. This software may be subject to other third-
@ -63,7 +63,7 @@ unsigned SpecAnalyzer::getMeanAbsValues (const int32_t* const mdctSignal, const
const unsigned channelIndex, const uint16_t* const bandStartOffsets, const unsigned nBands,
uint32_t* const meanBandValues)
{
if ((mdctSignal == nullptr) || (bandStartOffsets == nullptr) || (meanBandValues == nullptr) || (channelIndex >= USAC_MAX_NUM_CHANNELS) ||
if ((mdctSignal == nullptr) || (bandStartOffsets == nullptr) || (meanBandValues == nullptr) || (channelIndex > USAC_MAX_NUM_CHANNELS) ||
(nSamplesInFrame > 2048) || (nSamplesInFrame < 2) || (nBands > nSamplesInFrame))
{
return 1; // invalid arguments error
@ -77,7 +77,8 @@ unsigned SpecAnalyzer::getMeanAbsValues (const int32_t* const mdctSignal, const
const unsigned bandWidth = __min (nSamplesInFrame, bandStartOffsets[b + 1]) - bandOffset;
const unsigned anaBandIdx = bandOffset >> SA_BW_SHIFT;
if ((anaBandIdx < m_numAnaBands[channelIndex]) && (bandOffset == (anaBandIdx << SA_BW_SHIFT)) && ((bandWidth & (SA_BW - 1)) == 0))
if ((channelIndex < USAC_MAX_NUM_CHANNELS) && (anaBandIdx < m_numAnaBands[channelIndex]) &&
(bandOffset == (anaBandIdx << SA_BW_SHIFT)) && ((bandWidth & (SA_BW - 1)) == 0))
{
const uint32_t* const anaAbsVal = &m_meanAbsValue[channelIndex][anaBandIdx];
@ -94,13 +95,14 @@ unsigned SpecAnalyzer::getMeanAbsValues (const int32_t* const mdctSignal, const
{
#if SA_EXACT_COMPLEX_ABS
const double complexSqr = (double) bMdct[s] * (double) bMdct[s] + (double) bMdst[s] * (double) bMdst[s];
const unsigned absSample = unsigned (sqrt (complexSqr) + 0.5);
sumAbsVal += uint64_t (sqrt (complexSqr) + 0.5);
#else
const unsigned absReal = abs (bMdct[s]); // Richard Lyons, 1997; en.wikipedia.org/
const unsigned absImag = abs (bMdst[s]); // wiki/Alpha_max_plus_beta_min_algorithm
const unsigned absSample = (absReal > absImag ? absReal + ((absImag * 3) >> 3) : absImag + ((absReal * 3) >> 3));
const uint32_t absReal = abs (bMdct[s]); // Richard Lyons, 1997; en.wikipedia.org/
const uint32_t absImag = abs (bMdst[s]); // wiki/Alpha_max_plus_beta_min_algorithm
sumAbsVal += (absReal > absImag ? absReal + ((absImag * 3) >> 3) : absImag + ((absReal * 3) >> 3));
#endif
sumAbsVal += absSample;
}
// average spectral sample magnitude across current band
meanBandValues[b] = uint32_t ((sumAbsVal + (bandWidth >> 1)) / bandWidth);
@ -114,14 +116,34 @@ unsigned SpecAnalyzer::getMeanAbsValues (const int32_t* const mdctSignal, const
const unsigned bandOffset = __min (nSamplesInFrame, bandStartOffsets[b]);
const unsigned bandWidth = __min (nSamplesInFrame, bandStartOffsets[b + 1]) - bandOffset;
const int32_t* const bMdct = &mdctSignal[bandOffset];
#if SA_IMPROVED_REAL_ABS
uint64_t sumAbsVal = (bandStartOffsets[b + 1] == nSamplesInFrame ? abs (bMdct[bandWidth - 1]) : 0);
for (int s = bandWidth - (bandStartOffsets[b + 1] == nSamplesInFrame ? 2 : 1); s >= 0; s--)
{
// based on S. Merdjani, L. Daudet, "Estimation of Frequency from MDCT-Encoded Files,"
// DAFx-03, 2003, http://www.eecs.qmul.ac.uk/legacy/dafx03/proceedings/pdfs/dafx01.pdf
const uint32_t absReal = abs (bMdct[s]); // Richard Lyons, 1997; see also code above
const uint32_t absEstIm = abs (bMdct[s + 1] - bMdct[s - 1]) >> 1; // s - 1 may be -1!
sumAbsVal += (absReal > absEstIm ? absReal + ((absEstIm * 3) >> 3) : absEstIm + ((absReal * 3) >> 3));
}
if ((b == 0) && (bandWidth > 0)) // correct estimate at DC
{
const uint32_t absReal = abs (bMdct[0]); // Richard Lyons, 1997; see also code above
const uint32_t absEstIm = abs (bMdct[1] - bMdct[-1]) >> 1; // if s - 1 was -1 earlier
sumAbsVal -= (absReal > absEstIm ? absReal + ((absEstIm * 3) >> 3) : absEstIm + ((absReal * 3) >> 3));
sumAbsVal += absReal;
}
#else
uint64_t sumAbsVal = 0;
for (int s = bandWidth - 1; s >= 0; s--)
{
const unsigned absSample = abs (bMdct[s]);
sumAbsVal += absSample;
sumAbsVal += abs (bMdct[s]);
}
#endif
// average spectral sample magnitude across frequency band
meanBandValues[b] = uint32_t ((sumAbsVal + (bandWidth >> 1)) / bandWidth);
} // for b
@ -262,11 +284,11 @@ unsigned SpecAnalyzer::spectralAnalysis (const int32_t* const mdctSignals[USAC_M
// sum absolute values of complex signal, derive L1 norm, peak value, and peak index
#if SA_EXACT_COMPLEX_ABS
const double complexSqr = (double) bMdct[s] * (double) bMdct[s] + (double) bMdst[s] * (double) bMdst[s];
const unsigned absSample = unsigned (sqrt (complexSqr) + 0.5);
const uint32_t absSample = uint32_t (sqrt (complexSqr) + 0.5);
#else
const unsigned absReal = abs (bMdct[s]); // Richard Lyons, 1997; en.wikipedia.org/
const unsigned absImag = abs (bMdst[s]); // wiki/Alpha_max_plus_beta_min_algorithm
const unsigned absSample = (absReal > absImag ? absReal + ((absImag * 3) >> 3) : absImag + ((absReal * 3) >> 3));
const uint32_t absReal = abs (bMdct[s]); // Richard Lyons, 1997; en.wikipedia.org/
const uint32_t absImag = abs (bMdst[s]); // wiki/Alpha_max_plus_beta_min_algorithm
const uint32_t absSample = (absReal > absImag ? absReal + ((absImag * 3) >> 3) : absImag + ((absReal * 3) >> 3));
#endif
sumAbsVal += absSample;
if (offs + s > 0) // exclude DC from max/min
@ -288,7 +310,7 @@ unsigned SpecAnalyzer::spectralAnalysis (const int32_t* const mdctSignals[USAC_M
for (int s = SA_BW - 1; s >= 0; s--)
{
// obtain absolute values of real signal, derive L1 norm, peak value, and peak index
const unsigned absSample = abs (bMdct[s]);
const uint32_t absSample = abs (bMdct[s]);
sumAbsVal += absSample;
if (offs + s > 0) // exclude DC from max/min

View File

@ -1,5 +1,5 @@
/* specAnalysis.h - header file for class providing spectral analysis of MCLT signals
* written by C. R. Helmrich, last modified in 2019 - see License.htm for legal notices
* written by C. R. Helmrich, last modified in 2020 - see License.htm for legal notices
*
* The copyright in this software is being made available under a Modified BSD-Style License
* and comes with ABSOLUTELY NO WARRANTY. This software may be subject to other third-
@ -19,6 +19,7 @@
#define SA_BW (1 << SA_BW_SHIFT)
#define SA_EPS 1024
#define SA_EXACT_COMPLEX_ABS 0
#define SA_IMPROVED_REAL_ABS 1
#define SA_OPT_WINDOW_GROUPING 1
// spectral signal analysis class