From f12c25b7ae11bf4196292a6db3ed65bc668509a1 Mon Sep 17 00:00:00 2001 From: Max Kellermann Date: Fri, 3 Dec 2021 16:32:15 +0100 Subject: [PATCH] pcm/ReplayGainAnalyzer: new library --- src/pcm/ReplayGainAnalyzer.cxx | 363 +++++++++++++++++++++++++++++++++ src/pcm/ReplayGainAnalyzer.hxx | 142 +++++++++++++ src/pcm/meson.build | 1 + test/RunReplayGainAnalyzer.cxx | 63 ++++++ test/meson.build | 11 + 5 files changed, 580 insertions(+) create mode 100644 src/pcm/ReplayGainAnalyzer.cxx create mode 100644 src/pcm/ReplayGainAnalyzer.hxx create mode 100644 test/RunReplayGainAnalyzer.cxx diff --git a/src/pcm/ReplayGainAnalyzer.cxx b/src/pcm/ReplayGainAnalyzer.cxx new file mode 100644 index 000000000..10476f2c3 --- /dev/null +++ b/src/pcm/ReplayGainAnalyzer.cxx @@ -0,0 +1,363 @@ +/* + * Copyright 2021 Max Kellermann + * + * This library is based on af_replaygain.c from FFmpeg. Original + * copyright header: + * + * Copyright (c) 1998 - 2009 Conifer Software + * + * This file is part of FFmpeg. + * + * FFmpeg is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * FFmpeg is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with FFmpeg; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#include "ReplayGainAnalyzer.hxx" +#include "util/Compiler.h" +#include "util/ConstBuffer.hxx" + +#include +#include +#include +#include +#include + +ReplayGainAnalyzer::ReplayGainAnalyzer() noexcept +{ +} + +/* + * Find the largest absolute sample value. + */ +[[gnu::pure]] [[gnu::hot]] +static float +FindPeak(const float *samples, std::size_t n) noexcept +{ + float peak = 0.0; + + while (n-- > 0) { + float value = std::fabs(*samples++); + if (value > peak) + peak = value; + } + + return peak; +} + +static constexpr double +Square(double x) noexcept +{ + return x * x; +} + +[[gnu::const]] +static double +SquareHypot(ReplayGainAnalyzer::Frame f) noexcept +{ +#if GCC_OLDER_THAN(10,0) + /* GCC 8 doesn't have std::transform_reduce() */ + double sum = 0; + for (const auto &i : f) + sum += Square(i); + return sum; +#else + /* proper C++17 */ + return std::transform_reduce(f.begin(), f.end(), 0., + std::plus{}, Square); +#endif +} + +/* + * Calculate stereo RMS level. Minimum value is about -100 dB for + * digital silence. The 90 dB offset is to compensate for the + * normalized float range and 3 dB is for stereo samples. + */ +[[gnu::hot]] +static double +CalcStereoRMS(ConstBuffer src) noexcept +{ +#if GCC_OLDER_THAN(10,0) + /* GCC 8 doesn't have std::transform_reduce() */ + double sum = 0; + for (const auto &i : src) + sum += SquareHypot(i); +#else + /* proper C++17 */ + double sum = std::transform_reduce(src.begin(), src.end(), + 1e-16, + std::plus{}, + SquareHypot); +#endif + + return 10 * std::log10(sum / src.size) + 90.0 - 3.0; +} + +static constexpr bool +IsSilentSample(float value) noexcept +{ + return std::fabs(value) <= 1e-10f; +} + +[[gnu::const]] +static bool +IsSilentFrame(ReplayGainAnalyzer::Frame frame) noexcept +{ + return std::all_of(frame.begin(), frame.end(), IsSilentSample); +} + +[[gnu::pure]] +static bool +IsSilentBuffer(ConstBuffer buffer) noexcept +{ + return std::all_of(buffer.begin(), buffer.end(), IsSilentFrame); +} + +static constexpr ReplayGainAnalyzer::DoubleFrame +operator*(ReplayGainAnalyzer::Frame f, float x) noexcept +{ + ReplayGainAnalyzer::DoubleFrame result{}; + for (std::size_t i = 0; i < result.size(); ++i) + result[i] = (double)f[i] * (double)x; + return result; +} + +static constexpr auto & +operator+=(ReplayGainAnalyzer::DoubleFrame &dest, ReplayGainAnalyzer::DoubleFrame src) noexcept +{ + for (std::size_t i = 0; i < dest.size(); ++i) + dest[i] += src[i]; + return dest; +} + +static constexpr auto & +operator-=(ReplayGainAnalyzer::DoubleFrame &dest, ReplayGainAnalyzer::DoubleFrame src) noexcept +{ + for (std::size_t i = 0; i < dest.size(); ++i) + dest[i] -= src[i]; + return dest; +} + +[[maybe_unused]] +static constexpr auto +operator+(ReplayGainAnalyzer::DoubleFrame a, ReplayGainAnalyzer::DoubleFrame b) noexcept +{ + a += b; + return a; +} + +static constexpr auto +operator-(ReplayGainAnalyzer::DoubleFrame a, ReplayGainAnalyzer::DoubleFrame b) noexcept +{ + a -= b; + return a; +} + +static constexpr auto +ToSingle(ReplayGainAnalyzer::DoubleFrame src) noexcept +{ + ReplayGainAnalyzer::Frame dest{}; + for (std::size_t i = 0; i < dest.size(); ++i) + dest[i] = src[i]; + return dest; +} + +template +static constexpr ReplayGainAnalyzer::Frame +ApplyFilter(ReplayGainAnalyzer::Frame src, std::size_t i, + std::array &hist_a, + std::array &hist_b, + const std::array &coeff_a, + const std::array &coeff_b) noexcept +{ + ReplayGainAnalyzer::DoubleFrame frame = (hist_b[i] = src) * coeff_b[0]; + +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC unroll 100 +#endif + for (std::size_t j = 1; j <= ORDER; ++j) + frame += hist_b[i - j] * coeff_b[j] - hist_a[i - j] * coeff_a[j]; + + return hist_a[i] = ToSingle(frame); +} + +[[gnu::hot]] +inline void +ReplayGainAnalyzer::Yule::Filter(const Frame *gcc_restrict src, + Frame *gcc_restrict dst, + std::size_t n_frames) noexcept +{ + std::size_t i = hist_i; + + // If filter history is very small magnitude, clear it completely to + // prevent denormals from rattling around in there forever + // (slowing us down). + + if (IsSilentBuffer({&hist_a[i - ORDER], ORDER}) || + IsSilentBuffer({&hist_b[i - ORDER], ORDER})) + hist_a = hist_b = {}; + + while (n_frames--) { + *dst++ = ApplyFilter(*src++, i++, + hist_a, hist_b, + coeff_a, coeff_b); + + if (i == hist_a.size()) { + constexpr std::size_t n = ORDER; + std::copy(std::prev(hist_a.end(), n), hist_a.end(), + hist_a.begin()); + std::copy(std::prev(hist_b.end(), n), hist_b.end(), + hist_b.begin()); + i = n; + } + } + + hist_i = i; +} + +[[gnu::hot]] +inline void +ReplayGainAnalyzer::Butter::Filter(Frame *gcc_restrict samples, + std::size_t n_frames) noexcept +{ + std::size_t i = hist_i; + + // If filter history is very small magnitude, clear it completely + // to prevent denormals from rattling around in there forever + // (slowing us down). + + if (IsSilentBuffer({&hist_a[i - ORDER], ORDER}) || + IsSilentBuffer({&hist_b[i - ORDER], ORDER})) + hist_a = hist_b = {}; + + while (n_frames--) { + *samples = ApplyFilter(*samples, i, + hist_a, hist_b, + coeff_a, coeff_b); + ++samples; + ++i; + + if (i == hist_a.size()) { + constexpr std::size_t n = ORDER; + std::copy(std::prev(hist_a.end(), n), hist_a.end(), + hist_a.begin()); + std::copy(std::prev(hist_b.end(), n), hist_b.end(), + hist_b.begin()); + i = n; + } + } + + hist_i = i; +} + +void +ReplayGainAnalyzer::Process(ConstBuffer src) noexcept +{ + assert(!src.empty()); + + float new_peak = FindPeak(src.front().data(), + src.size * src.front().size()); + if (new_peak > peak) + peak = new_peak; + + Frame *tmp = buffer.GetT(src.size); + yule.Filter(src.data, tmp, src.size); + butter.Filter(tmp, src.size); + + const long level = std::lrint(std::floor(STEPS_PER_DB * CalcStereoRMS({tmp, src.size}))); + const std::size_t level_index = std::clamp(level, 0L, (long)histogram.size() - 1L); + histogram[level_index]++; +} + +/* + * Calculate the ReplayGain value from the specified loudness histogram; + * clip to -24 / +64 dB. + */ +template +[[gnu::pure]] [[gnu::hot]] +static float +FindHistogramPercentile95(const std::array &histogram) noexcept +{ +#if GCC_OLDER_THAN(10,0) + /* GCC 8 doesn't have std::reduce() */ + uint_fast32_t total_windows = 0; + for (const auto &i : histogram) + total_windows += i; +#else + /* proper C++17 */ + const uint_fast32_t total_windows = + std::reduce(histogram.begin(), histogram.end()); +#endif + + uint_fast32_t loud_count = 0; + std::size_t i = histogram.size(); + while (i--) + if ((loud_count += histogram[i]) * 20 >= total_windows) + break; + + return i; +} + +float +ReplayGainAnalyzer::GetGain() const noexcept +{ + std::size_t i = FindHistogramPercentile95(histogram); + float gain = (float)(64.54f - float(i) / STEPS_PER_DB); + + return std::clamp(gain, -24.0f, 64.0f); +} + +void +WindowReplayGainAnalyzer::CopyToBuffer(ConstBuffer src) noexcept +{ + std::copy(src.begin(), src.end(), + window_buffer.data() + window_fill); + window_fill += src.size; +} + +void +WindowReplayGainAnalyzer::Process(ConstBuffer src) noexcept +{ + assert(window_fill < WINDOW_FRAMES); + + if (window_fill > 0) { + std::size_t window_space = WINDOW_FRAMES - window_fill; + + if (src.size < window_space) { + CopyToBuffer(src); + return; + } + + CopyToBuffer({src.data, window_space}); + Flush(); + + src.skip_front(window_space); + if (src.empty()) + return; + } + + while (src.size >= WINDOW_FRAMES) { + ReplayGainAnalyzer::Process({src.data, WINDOW_FRAMES}); + src.skip_front(WINDOW_FRAMES); + } + + CopyToBuffer(src); +} + +void +WindowReplayGainAnalyzer::Flush() noexcept +{ + if (window_fill > 0) + ReplayGainAnalyzer::Process({window_buffer.data(), window_fill}); + window_fill = 0; +} diff --git a/src/pcm/ReplayGainAnalyzer.hxx b/src/pcm/ReplayGainAnalyzer.hxx new file mode 100644 index 000000000..fcba9ad4b --- /dev/null +++ b/src/pcm/ReplayGainAnalyzer.hxx @@ -0,0 +1,142 @@ +/* + * Copyright 2021 Max Kellermann + * + * This library is based on af_replaygain.c from FFmpeg. Original + * copyright header: + * + * Copyright (c) 1998 - 2009 Conifer Software + * + * This file is part of FFmpeg. + * + * FFmpeg is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * FFmpeg is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with FFmpeg; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#pragma once + +#include "Buffer.hxx" + +#include +#include + +template struct ConstBuffer; + +/** + * Analyze a 44.1 kHz / stereo / float32 audio stream and calculate + * peak and ReplayGain values. + */ +class ReplayGainAnalyzer { +public: + static constexpr unsigned CHANNELS = 2; + static constexpr unsigned SAMPLE_RATE = 44100; + using sample_type = float; + using Frame = std::array; + using DoubleFrame = std::array; + +private: + /* + * Optimized implementation of 10th-order IIR stereo filter. + */ + struct Yule { + static constexpr std::size_t ORDER = 10; + + static constexpr std::array coeff_a{ + 1.00000000000000, -3.47845948550071, 6.36317777566148, + -8.54751527471874, 9.47693607801280, -8.81498681370155, + 6.85401540936998, -4.39470996079559, 2.19611684890774, + -0.75104302451432, 0.13149317958808, + }; + + static constexpr std::array coeff_b{ + 0.05418656406430, -0.02911007808948, -0.00848709379851, + -0.00851165645469, -0.00834990904936, 0.02245293253339, + -0.02596338512915, 0.01624864962975, -0.00240879051584, + 0.00674613682247, -0.00187763777362 + }; + + static_assert(coeff_a.size() == ORDER + 1); + static_assert(coeff_b.size() == ORDER + 1); + + unsigned hist_i = ORDER; + std::array hist_a{}, hist_b{}; + + void Filter(const Frame *src, Frame *dst, + std::size_t n_frames) noexcept; + }; + + /* + * Optimized implementation of 2nd-order IIR stereo filter. + */ + struct Butter { + static constexpr std::size_t ORDER = 2; + + static constexpr std::array coeff_a{ + 1.00000000000000, -1.96977855582618, 0.97022847566350, + }; + + static constexpr std::array coeff_b{ + 0.98500175787242, -1.97000351574484, 0.98500175787242, + }; + + static_assert(coeff_a.size() == ORDER + 1); + static_assert(coeff_b.size() == ORDER + 1); + + unsigned hist_i = ORDER; + std::array hist_a{}, hist_b{}; + + void Filter(Frame *samples, std::size_t n_frames) noexcept; + }; + + static constexpr std::size_t STEPS_PER_DB = 100; + static constexpr unsigned MAX_DB = 120; + + std::array histogram{}; + float peak = 0; + + Yule yule; + Butter butter; + + PcmBuffer buffer; + +public: + ReplayGainAnalyzer() noexcept; + + void Process(ConstBuffer src) noexcept; + + float GetPeak() const noexcept { + return peak; + } + + [[gnu::pure]] + float GetGain() const noexcept; +}; + +/** + * A variant of #ReplayGainAnalyzer which automatically calls + * Process() with windows of 50ms. + */ +class WindowReplayGainAnalyzer : public ReplayGainAnalyzer { + static constexpr std::size_t WINDOW_FRAMES = SAMPLE_RATE / 20; + + std::array window_buffer; + std::size_t window_fill = 0; + +public: + void Process(ConstBuffer src) noexcept; + + void Flush() noexcept; + +private: + void CopyToBuffer(ConstBuffer src) noexcept; +}; diff --git a/src/pcm/meson.build b/src/pcm/meson.build index 74b605b5a..6fa0d6bc0 100644 --- a/src/pcm/meson.build +++ b/src/pcm/meson.build @@ -47,6 +47,7 @@ pcm_sources = [ 'FallbackResampler.cxx', 'ConfiguredResampler.cxx', 'AudioCompress/compress.c', + 'ReplayGainAnalyzer.cxx', ] libsamplerate_dep = dependency('samplerate', version: '>= 0.1.3', required: get_option('libsamplerate')) diff --git a/test/RunReplayGainAnalyzer.cxx b/test/RunReplayGainAnalyzer.cxx new file mode 100644 index 000000000..f179b17ed --- /dev/null +++ b/test/RunReplayGainAnalyzer.cxx @@ -0,0 +1,63 @@ +/* + * Copyright 2003-2021 The Music Player Daemon Project + * http://www.musicpd.org + * + * 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 + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#include "ReadFrames.hxx" +#include "pcm/ReplayGainAnalyzer.hxx" +#include "io/FileDescriptor.hxx" +#include "system/Error.hxx" +#include "util/ConstBuffer.hxx" +#include "util/PrintException.hxx" + +#include +#include + +#include + +int +main(int, char **) noexcept +try { + WindowReplayGainAnalyzer a; + + constexpr std::size_t frame_size = ReplayGainAnalyzer::CHANNELS * + sizeof(ReplayGainAnalyzer::sample_type); + + const FileDescriptor input_fd(STDIN_FILENO); + + while (true) { + std::array buffer; + + size_t nbytes = ReadFrames(input_fd, + buffer.data(), sizeof(buffer), + frame_size); + if (nbytes == 0) + break; + + a.Process({buffer.data(), nbytes / frame_size}); + } + + a.Flush(); + + printf("gain = %+.2f dB\n", (double)a.GetGain()); + printf("peak = %.6f\n", (double)a.GetPeak()); + + return EXIT_SUCCESS; +} catch (...) { + PrintException(std::current_exception()); + return EXIT_FAILURE; +} diff --git a/test/meson.build b/test/meson.build index 3b9671410..01bcfd0ba 100644 --- a/test/meson.build +++ b/test/meson.build @@ -556,6 +556,17 @@ executable( ], ) +executable( + 'RunReplayGainAnalyzer', + 'RunReplayGainAnalyzer.cxx', + 'ReadFrames.cxx', + include_directories: inc, + dependencies: [ + pcm_dep, + io_dep, + ], +) + # # Encoder #