openMSX
BlipBuffer.cc
Go to the documentation of this file.
1#include "BlipBuffer.hh"
2#include "cstd.hh"
3#include "Math.hh"
4#include "narrow.hh"
5#include "ranges.hh"
6#include "xrange.hh"
7#include <algorithm>
8#include <array>
9#include <cassert>
10#include <iostream>
11
12namespace openmsx {
13
14// Number of samples in a (pre-calculated) impulse-response wave-form.
15static constexpr int BLIP_IMPULSE_WIDTH = 16;
16
17// Derived constants
18static constexpr int BLIP_RES = 1 << BlipBuffer::BLIP_PHASE_BITS;
19
20
21// Precalculated impulse table.
22static constexpr auto impulses = [] {
23 constexpr int HALF_SIZE = BLIP_RES / 2 * (BLIP_IMPULSE_WIDTH - 1);
24 std::array<double, BLIP_RES + HALF_SIZE + BLIP_RES> fImpulse = {};
25 std::span<double, HALF_SIZE> out = subspan<HALF_SIZE>(fImpulse, BLIP_RES);
26 std::span<double, BLIP_RES> end = subspan<BLIP_RES >(fImpulse, BLIP_RES + HALF_SIZE);
27
28 // generate sinc, apply hamming window
29 double overSample = ((4.5 / (BLIP_IMPULSE_WIDTH - 1)) + 0.85);
30 double to_angle = Math::pi / (2.0 * overSample * BLIP_RES);
31 double to_fraction = Math::pi / (2 * (HALF_SIZE - 1));
32 for (auto i : xrange(HALF_SIZE)) {
33 double angle = ((i - HALF_SIZE) * 2 + 1) * to_angle;
34 out[i] = cstd::sin<2>(angle) / angle;
35 out[i] *= 0.54 - 0.46 * cstd::cos<2>((2 * i + 1) * to_fraction);
36 }
37
38 // need mirror slightly past center for calculation
39 for (auto i : xrange(BLIP_RES)) {
40 end[i] = out[HALF_SIZE - 1 - i];
41 }
42
43 // find rescale factor
44 double total = 0.0;
45 for (auto i : xrange(HALF_SIZE)) {
46 total += out[i];
47 }
48 double rescale = 1.0 / (2.0 * total);
49
50 // integrate, first difference, rescale, convert to float
51 constexpr int IMPULSES_SIZE = BLIP_RES * (BLIP_IMPULSE_WIDTH / 2) + 1;
52 std::array<float, IMPULSES_SIZE> imp = {};
53 double sum = 0.0;
54 double next = 0.0;
55 for (auto i : xrange(IMPULSES_SIZE)) {
56 imp[i] = float((next - sum) * rescale);
57 sum += fImpulse[i];
58 next += fImpulse[i + BLIP_RES];
59 }
60 // Original code would now apply a correction on each kernel so that
61 // the (integer) coefficients sum up to 'kernelUnit'. I've measured
62 // that after switching to float coefficients this correction is
63 // roughly equal to std::numeric_limits<float>::epsilon(). So it's no
64 // longer meaningful.
65
66 // reshuffle values to a more cache friendly order
67 std::array<std::array<float, BLIP_IMPULSE_WIDTH>, BLIP_RES> result = {};
68 for (auto phase : xrange(BLIP_RES)) {
69 const auto* imp_fwd = &imp[BLIP_RES - phase];
70 const auto* imp_rev = &imp[phase];
71 auto* p = result[phase].data();
72 for (size_t i : xrange(BLIP_IMPULSE_WIDTH / 2)) {
73 *p++ = imp_fwd[BLIP_RES * i];
74 }
75 for (ptrdiff_t i = BLIP_IMPULSE_WIDTH / 2 - 1; i >= 0; --i) {
76 *p++ = imp_rev[BLIP_RES * i];
77 }
78 }
79 return result;
80}();
81
82
84{
85 if (false) {
86 for (auto i : xrange(BLIP_RES)) {
87 std::cout << "\t{ " << impulses[i][0];
88 for (auto j : xrange(1, BLIP_IMPULSE_WIDTH)) {
89 std::cout << ", " << impulses[i][j];
90 }
91 std::cout << " },\n";
92 }
93 }
94
95 ranges::fill(buffer, 0);
96}
97
98void BlipBuffer::addDelta(TimeIndex time, float delta)
99{
100 unsigned tmp = time.toInt() + BLIP_IMPULSE_WIDTH;
101 assert(tmp < BUFFER_SIZE);
102 availSamp = std::max(availSamp, narrow<ptrdiff_t>(tmp));
103
104 auto phase = time.fractAsInt();
105 auto ofst = time.toInt() + offset;
106 const float* __restrict impulse = impulses[phase].data();
107 if ((ofst + BLIP_IMPULSE_WIDTH) <= BUFFER_SIZE) [[likely]] {
108 float* __restrict result = &buffer[ofst];
109 for (auto i : xrange(BLIP_IMPULSE_WIDTH)) {
110 result[i] += impulse[i] * delta;
111 }
112 } else {
113 for (auto i : xrange(BLIP_IMPULSE_WIDTH)) {
114 buffer[(ofst + i) & BUFFER_MASK] += impulse[i] * delta;
115 }
116 }
117}
118
119static constexpr float BASS_FACTOR = 511.0f / 512.0f;
120
121template<size_t PITCH>
122void BlipBuffer::readSamplesHelper(float* __restrict out, size_t samples)
123{
124 assert((offset + samples) <= BUFFER_SIZE);
125 auto acc = accum;
126 auto ofst = offset;
127 for (auto i : xrange(samples)) {
128 out[i * PITCH] = acc;
129 acc *= BASS_FACTOR;
130 acc += buffer[ofst];
131 buffer[ofst] = 0.0f;
132 ++ofst;
133 }
134 accum = acc;
135 offset = ofst & BUFFER_MASK;
136}
137
138static bool isSilent(float x)
139{
140 // 'accum' falls away slowly (via BASS_FACTOR). Because it's a float it
141 // takes a long time before it's really zero. But much sooner we can
142 // already say it's practically silent. This check is safe when the
143 // input is in range [-1..+1] (does not say silent too soon). When the
144 // input is in range [-32768..+32767] it takes longer before we switch
145 // to silent mode (but still less than a second).
146 constexpr float threshold = 1.0f / 32768;
147 return std::abs(x) < threshold;
148}
149
150// TODO replace 'out + samples + PITCH' with 'stride_view' ???
151template<size_t PITCH>
152bool BlipBuffer::readSamples(float* __restrict out, size_t samples)
153{
154 if (availSamp <= 0) {
155 #ifdef DEBUG
156 // buffer contains all zeros (only check this in debug mode)
157 assert(ranges::all_of(buffer, [](const auto& b) { return b == 0.0f; }));
158 #endif
159 if (isSilent(accum)) {
160 return false; // muted
161 }
162 auto acc = accum;
163 for (auto i : xrange(samples)) {
164 out[i * PITCH] = acc;
165 acc *= BASS_FACTOR;
166 }
167 accum = acc;
168 } else {
169 availSamp -= narrow<ptrdiff_t>(samples);
170 auto t1 = std::min(samples, BUFFER_SIZE - offset);
171 readSamplesHelper<PITCH>(out, t1);
172 if (t1 < samples) {
173 assert(offset == 0);
174 auto t2 = samples - t1;
175 assert(t2 < BUFFER_SIZE);
176 readSamplesHelper<PITCH>(&out[t1 * PITCH], t2);
177 }
178 assert(offset < BUFFER_SIZE);
179 }
180 return true;
181}
182
183template bool BlipBuffer::readSamples<1>(float*, size_t);
184template bool BlipBuffer::readSamples<2>(float*, size_t);
185
186} // namespace openmsx
bool readSamples(float *out, size_t samples)
Definition: BlipBuffer.cc:152
static constexpr int BLIP_PHASE_BITS
Definition: BlipBuffer.hh:20
void addDelta(TimeIndex time, float delta)
Definition: BlipBuffer.cc:98
A fixed point number, implemented by a 32-bit signed integer.
Definition: FixedPoint.hh:15
constexpr int toInt() const
Returns the integer part (rounded down) of this fixed point number.
Definition: FixedPoint.hh:76
constexpr unsigned fractAsInt() const
Returns the fractional part of this value as an integer.
Definition: FixedPoint.hh:144
constexpr double pi
Definition: Math.hh:23
constexpr vecN< N, T > min(const vecN< N, T > &x, const vecN< N, T > &y)
Definition: gl_vec.hh:266
constexpr vecN< N, T > max(const vecN< N, T > &x, const vecN< N, T > &y)
Definition: gl_vec.hh:284
This file implemented 3 utility functions:
Definition: Autofire.cc:9
bool all_of(InputRange &&range, UnaryPredicate pred)
Definition: ranges.hh:186
constexpr void fill(ForwardRange &&range, const T &value)
Definition: ranges.hh:287
uint32_t next(octet_iterator &it, octet_iterator end)
auto sum(InputRange &&range, Proj proj={})
Definition: stl.hh:236
constexpr auto xrange(T e)
Definition: xrange.hh:133
constexpr auto end(const zstring_view &x)