| 1 | #include "./common.h" |
| 2 | |
| 3 | #ifndef SIGNALSMITH_DSP_DELAY_H |
| 4 | #define SIGNALSMITH_DSP_DELAY_H |
| 5 | |
| 6 | #include <vector> |
| 7 | #include <array> |
| 8 | #include <cmath> // for std::ceil() |
| 9 | #include <type_traits> |
| 10 | |
| 11 | #include <complex> |
| 12 | #include "./fft.h" |
| 13 | #include "./windows.h" |
| 14 | |
| 15 | namespace signalsmith { |
| 16 | namespace delay { |
| 17 | /** @defgroup Delay Delay utilities |
| 18 | @brief Standalone templated classes for delays |
| 19 | |
| 20 | You can set up a `Buffer` or `MultiBuffer`, and get interpolated samples using a `Reader` (separately on each channel in the multi-channel case) - or you can use `Delay`/`MultiDelay` which include their own buffers. |
| 21 | |
| 22 | Interpolation quality is chosen using a template class, from @ref Interpolators. |
| 23 | |
| 24 | @{ |
| 25 | @file |
| 26 | */ |
| 27 | |
| 28 | /** @brief Single-channel delay buffer |
| 29 | |
| 30 | Access is used with `buffer[]`, relative to the internal read/write position ("head"). This head is moved using `++buffer` (or `buffer += n`), such that `buffer[1] == (buffer + 1)[0]` in a similar way iterators/pointers. |
| 31 | |
| 32 | Operations like `buffer - 10` or `buffer++` return a View, which holds a fixed position in the buffer (based on the read/write position at the time). |
| 33 | |
| 34 | The capacity includes both positive and negative indices. For example, a capacity of 100 would support using any of the ranges: |
| 35 | |
| 36 | * `buffer[-99]` to buffer[0]` |
| 37 | * `buffer[-50]` to buffer[49]` |
| 38 | * `buffer[0]` to buffer[99]` |
| 39 | |
| 40 | Although buffers are usually used with historical samples accessed using negative indices e.g. `buffer[-10]`, you could equally use it flipped around (moving the head backwards through the buffer using `--buffer`). |
| 41 | */ |
| 42 | template<typename Sample> |
| 43 | class Buffer { |
| 44 | unsigned bufferIndex; |
| 45 | unsigned bufferMask; |
| 46 | std::vector<Sample> buffer; |
| 47 | public: |
| 48 | Buffer(int minCapacity=0) { |
| 49 | resize(minCapacity); |
| 50 | } |
| 51 | // We shouldn't accidentally copy a delay buffer |
| 52 | Buffer(const Buffer &other) = delete; |
| 53 | Buffer & operator =(const Buffer &other) = delete; |
| 54 | // But moving one is fine |
| 55 | Buffer(Buffer &&other) = default; |
| 56 | Buffer & operator =(Buffer &&other) = default; |
| 57 | |
| 58 | void resize(int minCapacity, Sample value=Sample()) { |
| 59 | int bufferLength = 1; |
| 60 | while (bufferLength < minCapacity) bufferLength *= 2; |
| 61 | buffer.assign(bufferLength, value); |
| 62 | bufferMask = unsigned(bufferLength - 1); |
| 63 | bufferIndex = 0; |
| 64 | } |
| 65 | void reset(Sample value=Sample()) { |
| 66 | buffer.assign(buffer.size(), value); |
| 67 | } |
| 68 | |
| 69 | /// Holds a view for a particular position in the buffer |
| 70 | template<bool isConst> |
| 71 | class View { |
| 72 | using CBuffer = typename std::conditional<isConst, const Buffer, Buffer>::type; |
| 73 | using CSample = typename std::conditional<isConst, const Sample, Sample>::type; |
| 74 | CBuffer *buffer = nullptr; |
| 75 | unsigned bufferIndex = 0; |
| 76 | public: |
| 77 | View(CBuffer &buffer, int offset=0) : buffer(&buffer), bufferIndex(buffer.bufferIndex + (unsigned)offset) {} |
| 78 | View(const View &other, int offset=0) : buffer(other.buffer), bufferIndex(other.bufferIndex + (unsigned)offset) {} |
| 79 | View & operator =(const View &other) { |
| 80 | buffer = other.buffer; |
| 81 | bufferIndex = other.bufferIndex; |
| 82 | return *this; |
| 83 | } |
| 84 | |
| 85 | CSample & operator[](int offset) { |
| 86 | return buffer->buffer[(bufferIndex + (unsigned)offset)&buffer->bufferMask]; |
| 87 | } |
| 88 | const Sample & operator[](int offset) const { |
| 89 | return buffer->buffer[(bufferIndex + (unsigned)offset)&buffer->bufferMask]; |
| 90 | } |
| 91 | |
| 92 | /// Write data into the buffer |
| 93 | template<typename Data> |
| 94 | void write(Data &&data, int length) { |
| 95 | for (int i = 0; i < length; ++i) { |
| 96 | (*this)[i] = data[i]; |
| 97 | } |
| 98 | } |
| 99 | /// Read data out from the buffer |
| 100 | template<typename Data> |
| 101 | void read(int length, Data &&data) const { |
| 102 | for (int i = 0; i < length; ++i) { |
| 103 | data[i] = (*this)[i]; |
| 104 | } |
| 105 | } |
| 106 | |
| 107 | View operator +(int offset) const { |
| 108 | return View(*this, offset); |
| 109 | } |
| 110 | View operator -(int offset) const { |
| 111 | return View(*this, -offset); |
| 112 | } |
| 113 | }; |
| 114 | using MutableView = View<false>; |
| 115 | using ConstView = View<true>; |
| 116 | |
| 117 | MutableView view(int offset=0) { |
| 118 | return MutableView(*this, offset); |
| 119 | } |
| 120 | ConstView view(int offset=0) const { |
| 121 | return ConstView(*this, offset); |
| 122 | } |
| 123 | ConstView constView(int offset=0) const { |
| 124 | return ConstView(*this, offset); |
| 125 | } |
| 126 | |
| 127 | Sample & operator[](int offset) { |
| 128 | return buffer[(bufferIndex + (unsigned)offset)&bufferMask]; |
| 129 | } |
| 130 | const Sample & operator[](int offset) const { |
| 131 | return buffer[(bufferIndex + (unsigned)offset)&bufferMask]; |
| 132 | } |
| 133 | |
| 134 | /// Write data into the buffer |
| 135 | template<typename Data> |
| 136 | void write(Data &&data, int length) { |
| 137 | for (int i = 0; i < length; ++i) { |
| 138 | (*this)[i] = data[i]; |
| 139 | } |
| 140 | } |
| 141 | /// Read data out from the buffer |
| 142 | template<typename Data> |
| 143 | void read(int length, Data &&data) const { |
| 144 | for (int i = 0; i < length; ++i) { |
| 145 | data[i] = (*this)[i]; |
| 146 | } |
| 147 | } |
| 148 | |
| 149 | Buffer & operator ++() { |
| 150 | ++bufferIndex; |
| 151 | return *this; |
| 152 | } |
| 153 | Buffer & operator +=(int i) { |
| 154 | bufferIndex += (unsigned)i; |
| 155 | return *this; |
| 156 | } |
| 157 | Buffer & operator --() { |
| 158 | --bufferIndex; |
| 159 | return *this; |
| 160 | } |
| 161 | Buffer & operator -=(int i) { |
| 162 | bufferIndex -= (unsigned)i; |
| 163 | return *this; |
| 164 | } |
| 165 | |
| 166 | MutableView operator ++(int) { |
| 167 | MutableView view(*this); |
| 168 | ++bufferIndex; |
| 169 | return view; |
| 170 | } |
| 171 | MutableView operator +(int i) { |
| 172 | return MutableView(*this, i); |
| 173 | } |
| 174 | ConstView operator +(int i) const { |
| 175 | return ConstView(*this, i); |
| 176 | } |
| 177 | MutableView operator --(int) { |
| 178 | MutableView view(*this); |
| 179 | --bufferIndex; |
| 180 | return view; |
| 181 | } |
| 182 | MutableView operator -(int i) { |
| 183 | return MutableView(*this, -i); |
| 184 | } |
| 185 | ConstView operator -(int i) const { |
| 186 | return ConstView(*this, -i); |
| 187 | } |
| 188 | }; |
| 189 | |
| 190 | /** @brief Multi-channel delay buffer |
| 191 | |
| 192 | This behaves similarly to the single-channel `Buffer`, with the following differences: |
| 193 | |
| 194 | * `buffer[c]` returns a view for a single channel, which behaves like the single-channel `Buffer::View`. |
| 195 | * The constructor and `.resize()` take an additional first `channel` argument. |
| 196 | */ |
| 197 | template<typename Sample> |
| 198 | class MultiBuffer { |
| 199 | int channels, stride; |
| 200 | Buffer<Sample> buffer; |
| 201 | public: |
| 202 | using ConstChannel = typename Buffer<Sample>::ConstView; |
| 203 | using MutableChannel = typename Buffer<Sample>::MutableView; |
| 204 | |
| 205 | MultiBuffer(int channels=0, int capacity=0) : channels(channels), stride(capacity), buffer(channels*capacity) {} |
| 206 | |
| 207 | void resize(int nChannels, int capacity, Sample value=Sample()) { |
| 208 | channels = nChannels; |
| 209 | stride = capacity; |
| 210 | buffer.resize(channels*capacity, value); |
| 211 | } |
| 212 | void reset(Sample value=Sample()) { |
| 213 | buffer.reset(value); |
| 214 | } |
| 215 | |
| 216 | /// A reference-like multi-channel result for a particular sample index |
| 217 | template<bool isConst> |
| 218 | class Stride { |
| 219 | using CChannel = typename std::conditional<isConst, ConstChannel, MutableChannel>::type; |
| 220 | using CSample = typename std::conditional<isConst, const Sample, Sample>::type; |
| 221 | CChannel view; |
| 222 | int channels, stride; |
| 223 | public: |
| 224 | Stride(CChannel view, int channels, int stride) : view(view), channels(channels), stride(stride) {} |
| 225 | Stride(const Stride &other) : view(other.view), channels(other.channels), stride(other.stride) {} |
| 226 | |
| 227 | CSample & operator[](int channel) { |
| 228 | return view[channel*stride]; |
| 229 | } |
| 230 | const Sample & operator[](int channel) const { |
| 231 | return view[channel*stride]; |
| 232 | } |
| 233 | |
| 234 | /// Reads from the buffer into a multi-channel result |
| 235 | template<class Data> |
| 236 | void get(Data &&result) const { |
| 237 | for (int c = 0; c < channels; ++c) { |
| 238 | result[c] = view[c*stride]; |
| 239 | } |
| 240 | } |
| 241 | /// Writes from multi-channel data into the buffer |
| 242 | template<class Data> |
| 243 | void set(Data &&data) { |
| 244 | for (int c = 0; c < channels; ++c) { |
| 245 | view[c*stride] = data[c]; |
| 246 | } |
| 247 | } |
| 248 | template<class Data> |
| 249 | Stride & operator =(const Data &data) { |
| 250 | set(data); |
| 251 | return *this; |
| 252 | } |
| 253 | Stride & operator =(const Stride &data) { |
| 254 | set(data); |
| 255 | return *this; |
| 256 | } |
| 257 | }; |
| 258 | |
| 259 | Stride<false> at(int offset) { |
| 260 | return {buffer.view(offset), channels, stride}; |
| 261 | } |
| 262 | Stride<true> at(int offset) const { |
| 263 | return {buffer.view(offset), channels, stride}; |
| 264 | } |
| 265 | |
| 266 | /// Holds a particular position in the buffer |
| 267 | template<bool isConst> |
| 268 | class View { |
| 269 | using CChannel = typename std::conditional<isConst, ConstChannel, MutableChannel>::type; |
| 270 | CChannel view; |
| 271 | int channels, stride; |
| 272 | public: |
| 273 | View(CChannel view, int channels, int stride) : view(view), channels(channels), stride(stride) {} |
| 274 | |
| 275 | CChannel operator[](int channel) { |
| 276 | return view + channel*stride; |
| 277 | } |
| 278 | ConstChannel operator[](int channel) const { |
| 279 | return view + channel*stride; |
| 280 | } |
| 281 | |
| 282 | Stride<isConst> at(int offset) { |
| 283 | return {view + offset, channels, stride}; |
| 284 | } |
| 285 | Stride<true> at(int offset) const { |
| 286 | return {view + offset, channels, stride}; |
| 287 | } |
| 288 | }; |
| 289 | using MutableView = View<false>; |
| 290 | using ConstView = View<true>; |
| 291 | |
| 292 | MutableView view(int offset=0) { |
| 293 | return MutableView(buffer.view(offset), channels, stride); |
| 294 | } |
| 295 | ConstView view(int offset=0) const { |
| 296 | return ConstView(buffer.view(offset), channels, stride); |
| 297 | } |
| 298 | ConstView constView(int offset=0) const { |
| 299 | return ConstView(buffer.view(offset), channels, stride); |
| 300 | } |
| 301 | |
| 302 | MutableChannel operator[](int channel) { |
| 303 | return buffer + channel*stride; |
| 304 | } |
| 305 | ConstChannel operator[](int channel) const { |
| 306 | return buffer + channel*stride; |
| 307 | } |
| 308 | |
| 309 | MultiBuffer & operator ++() { |
| 310 | ++buffer; |
| 311 | return *this; |
| 312 | } |
| 313 | MultiBuffer & operator +=(int i) { |
| 314 | buffer += i; |
| 315 | return *this; |
| 316 | } |
| 317 | MutableView operator ++(int) { |
| 318 | return MutableView(buffer++, channels, stride); |
| 319 | } |
| 320 | MutableView operator +(int i) { |
| 321 | return MutableView(buffer + i, channels, stride); |
| 322 | } |
| 323 | ConstView operator +(int i) const { |
| 324 | return ConstView(buffer + i, channels, stride); |
| 325 | } |
| 326 | MultiBuffer & operator --() { |
| 327 | --buffer; |
| 328 | return *this; |
| 329 | } |
| 330 | MultiBuffer & operator -=(int i) { |
| 331 | buffer -= i; |
| 332 | return *this; |
| 333 | } |
| 334 | MutableView operator --(int) { |
| 335 | return MutableView(buffer--, channels, stride); |
| 336 | } |
| 337 | MutableView operator -(int i) { |
| 338 | return MutableView(buffer - i, channels, stride); |
| 339 | } |
| 340 | ConstView operator -(int i) const { |
| 341 | return ConstView(buffer - i, channels, stride); |
| 342 | } |
| 343 | }; |
| 344 | |
| 345 | /** \defgroup Interpolators Interpolators |
| 346 | \ingroup Delay |
| 347 | @{ */ |
| 348 | /// Nearest-neighbour interpolator |
| 349 | /// \diagram{delay-random-access-nearest.svg,aliasing and maximum amplitude/delay errors for different input frequencies} |
| 350 | template<typename Sample> |
| 351 | struct InterpolatorNearest { |
| 352 | static constexpr int inputLength = 1; |
| 353 | static constexpr Sample latency = -0.5; // Because we're truncating, which rounds down too often |
| 354 | |
| 355 | template<class Data> |
| 356 | static Sample fractional(const Data &data, Sample) { |
| 357 | return data[0]; |
| 358 | } |
| 359 | }; |
| 360 | /// Linear interpolator |
| 361 | /// \diagram{delay-random-access-linear.svg,aliasing and maximum amplitude/delay errors for different input frequencies} |
| 362 | template<typename Sample> |
| 363 | struct InterpolatorLinear { |
| 364 | static constexpr int inputLength = 2; |
| 365 | static constexpr int latency = 0; |
| 366 | |
| 367 | template<class Data> |
| 368 | static Sample fractional(const Data &data, Sample fractional) { |
| 369 | Sample a = data[0], b = data[1]; |
| 370 | return a + fractional*(b - a); |
| 371 | } |
| 372 | }; |
| 373 | /// Spline cubic interpolator |
| 374 | /// \diagram{delay-random-access-cubic.svg,aliasing and maximum amplitude/delay errors for different input frequencies} |
| 375 | template<typename Sample> |
| 376 | struct InterpolatorCubic { |
| 377 | static constexpr int inputLength = 4; |
| 378 | static constexpr int latency = 1; |
| 379 | |
| 380 | template<class Data> |
| 381 | static Sample fractional(const Data &data, Sample fractional) { |
| 382 | // Cubic interpolation |
| 383 | Sample a = data[0], b = data[1], c = data[2], d = data[3]; |
| 384 | Sample cbDiff = c - b; |
| 385 | Sample k1 = (c - a)*0.5; |
| 386 | Sample k3 = k1 + (d - b)*0.5 - cbDiff*2; |
| 387 | Sample k2 = cbDiff - k3 - k1; |
| 388 | return b + fractional*(k1 + fractional*(k2 + fractional*k3)); // 16 ops total, not including the indexing |
| 389 | } |
| 390 | }; |
| 391 | |
| 392 | // Efficient Algorithms and Structures for Fractional Delay Filtering Based on Lagrange Interpolation |
| 393 | // Franck 2009 https://www.aes.org/e-lib/browse.cfm?elib=14647 |
| 394 | namespace _franck_impl { |
| 395 | template<typename Sample, int n, int low, int high> |
| 396 | struct ProductRange { |
| 397 | using Array = std::array<Sample, (n + 1)>; |
| 398 | static constexpr int mid = (low + high)/2; |
| 399 | using Left = ProductRange<Sample, n, low, mid>; |
| 400 | using Right = ProductRange<Sample, n, mid + 1, high>; |
| 401 | |
| 402 | Left left; |
| 403 | Right right; |
| 404 | |
| 405 | const Sample total; |
| 406 | ProductRange(Sample x) : left(x), right(x), total(left.total*right.total) {} |
| 407 | |
| 408 | template<class Data> |
| 409 | Sample calculateResult(Sample , const Data &data, const Array &invFactors) { |
| 410 | return left.calculateResult(extraFactor*right.total, data, invFactors) |
| 411 | + right.calculateResult(extraFactor*left.total, data, invFactors); |
| 412 | } |
| 413 | }; |
| 414 | template<typename Sample, int n, int index> |
| 415 | struct ProductRange<Sample, n, index, index> { |
| 416 | using Array = std::array<Sample, (n + 1)>; |
| 417 | |
| 418 | const Sample total; |
| 419 | ProductRange(Sample x) : total(x - index) {} |
| 420 | |
| 421 | template<class Data> |
| 422 | Sample calculateResult(Sample , const Data &data, const Array &invFactors) { |
| 423 | return extraFactor*data[index]*invFactors[index]; |
| 424 | } |
| 425 | }; |
| 426 | } |
| 427 | /** Fixed-order Lagrange interpolation. |
| 428 | \diagram{interpolator-LagrangeN.svg,aliasing and amplitude/delay errors for different sizes} |
| 429 | */ |
| 430 | template<typename Sample, int n> |
| 431 | struct InterpolatorLagrangeN { |
| 432 | static constexpr int inputLength = n + 1; |
| 433 | static constexpr int latency = (n - 1)/2; |
| 434 | |
| 435 | using Array = std::array<Sample, (n + 1)>; |
| 436 | Array invDivisors; |
| 437 | |
| 438 | InterpolatorLagrangeN() { |
| 439 | for (int j = 0; j <= n; ++j) { |
| 440 | double divisor = 1; |
| 441 | for (int k = 0; k < j; ++k) divisor *= (j - k); |
| 442 | for (int k = j + 1; k <= n; ++k) divisor *= (j - k); |
| 443 | invDivisors[j] = 1/divisor; |
| 444 | } |
| 445 | } |
| 446 | |
| 447 | template<class Data> |
| 448 | Sample fractional(const Data &data, Sample fractional) const { |
| 449 | constexpr int mid = n/2; |
| 450 | using Left = _franck_impl::ProductRange<Sample, n, 0, mid>; |
| 451 | using Right = _franck_impl::ProductRange<Sample, n, mid + 1, n>; |
| 452 | |
| 453 | Sample x = fractional + latency; |
| 454 | |
| 455 | Left left(x); |
| 456 | Right right(x); |
| 457 | |
| 458 | return left.calculateResult(right.total, data, invDivisors) + right.calculateResult(left.total, data, invDivisors); |
| 459 | } |
| 460 | }; |
| 461 | template<typename Sample> |
| 462 | using InterpolatorLagrange3 = InterpolatorLagrangeN<Sample, 3>; |
| 463 | template<typename Sample> |
| 464 | using InterpolatorLagrange7 = InterpolatorLagrangeN<Sample, 7>; |
| 465 | template<typename Sample> |
| 466 | using InterpolatorLagrange19 = InterpolatorLagrangeN<Sample, 19>; |
| 467 | |
| 468 | /** Fixed-size Kaiser-windowed sinc interpolation. |
| 469 | \diagram{interpolator-KaiserSincN.svg,aliasing and amplitude/delay errors for different sizes} |
| 470 | If `minimumPhase` is enabled, a minimum-phase version of the kernel is used: |
| 471 | \diagram{interpolator-KaiserSincN-min.svg,aliasing and amplitude/delay errors for minimum-phase mode} |
| 472 | */ |
| 473 | template<typename Sample, int n, bool minimumPhase=false> |
| 474 | struct InterpolatorKaiserSincN { |
| 475 | static constexpr int inputLength = n; |
| 476 | static constexpr Sample latency = minimumPhase ? 0 : (n*Sample(0.5) - 1); |
| 477 | |
| 478 | int subSampleSteps; |
| 479 | std::vector<Sample> coefficients; |
| 480 | |
| 481 | InterpolatorKaiserSincN() : InterpolatorKaiserSincN(0.5 - 0.45/std::sqrt(x: n)) {} |
| 482 | InterpolatorKaiserSincN(double passFreq) : InterpolatorKaiserSincN(passFreq, 1 - passFreq) {} |
| 483 | InterpolatorKaiserSincN(double passFreq, double stopFreq) { |
| 484 | subSampleSteps = 2*n; // Heuristic again. Really it depends on the bandwidth as well. |
| 485 | double kaiserBandwidth = (stopFreq - passFreq)*(n + 1.0/subSampleSteps); |
| 486 | kaiserBandwidth += 1.25/kaiserBandwidth; // We want to place the first zero, but (because using this to window a sinc essentially integrates it in the freq-domain), our ripples (and therefore zeroes) are out of phase. This is a heuristic fix. |
| 487 | double sincScale = M_PI*(passFreq + stopFreq); |
| 488 | |
| 489 | double centreIndex = n*subSampleSteps*0.5, scaleFactor = 1.0/subSampleSteps; |
| 490 | std::vector<Sample> windowedSinc(subSampleSteps*n + 1); |
| 491 | |
| 492 | ::signalsmith::windows::Kaiser::withBandwidth(bandwidth: kaiserBandwidth, heuristicOptimal: false).fill(windowedSinc, windowedSinc.size()); |
| 493 | |
| 494 | for (size_t i = 0; i < windowedSinc.size(); ++i) { |
| 495 | double x = (i - centreIndex)*scaleFactor; |
| 496 | int intX = std::round(x: x); |
| 497 | if (intX != 0 && std::abs(x: x - intX) < 1e-6) { |
| 498 | // Exact 0s |
| 499 | windowedSinc[i] = 0; |
| 500 | } else if (std::abs(x: x) > 1e-6) { |
| 501 | double p = x*sincScale; |
| 502 | windowedSinc[i] *= std::sin(x: p)/p; |
| 503 | } |
| 504 | } |
| 505 | |
| 506 | if (minimumPhase) { |
| 507 | signalsmith::fft::FFT<Sample> fft(windowedSinc.size()*2, 1); |
| 508 | windowedSinc.resize(fft.size(), 0); |
| 509 | std::vector<std::complex<Sample>> spectrum(fft.size()); |
| 510 | std::vector<std::complex<Sample>> cepstrum(fft.size()); |
| 511 | fft.fft(windowedSinc, spectrum); |
| 512 | for (size_t i = 0; i < fft.size(); ++i) { |
| 513 | spectrum[i] = std::log(std::abs(spectrum[i]) + 1e-30); |
| 514 | } |
| 515 | fft.fft(spectrum, cepstrum); |
| 516 | for (size_t i = 1; i < fft.size()/2; ++i) { |
| 517 | cepstrum[i] *= 0; |
| 518 | } |
| 519 | for (size_t i = fft.size()/2 + 1; i < fft.size(); ++i) { |
| 520 | cepstrum[i] *= 2; |
| 521 | } |
| 522 | Sample scaling = Sample(1)/fft.size(); |
| 523 | fft.ifft(cepstrum, spectrum); |
| 524 | |
| 525 | for (size_t i = 0; i < fft.size(); ++i) { |
| 526 | Sample phase = spectrum[i].imag()*scaling; |
| 527 | Sample mag = std::exp(spectrum[i].real()*scaling); |
| 528 | spectrum[i] = {mag*std::cos(phase), mag*std::sin(phase)}; |
| 529 | } |
| 530 | fft.ifft(spectrum, cepstrum); |
| 531 | windowedSinc.resize(subSampleSteps*n + 1); |
| 532 | windowedSinc.shrink_to_fit(); |
| 533 | for (size_t i = 0; i < windowedSinc.size(); ++i) { |
| 534 | windowedSinc[i] = cepstrum[i].real()*scaling; |
| 535 | } |
| 536 | } |
| 537 | |
| 538 | // Re-order into FIR fractional-delay blocks |
| 539 | coefficients.resize(n*(subSampleSteps + 1)); |
| 540 | for (int k = 0; k <= subSampleSteps; ++k) { |
| 541 | for (int i = 0; i < n; ++i) { |
| 542 | coefficients[k*n + i] = windowedSinc[(subSampleSteps - k) + i*subSampleSteps]; |
| 543 | } |
| 544 | } |
| 545 | } |
| 546 | |
| 547 | template<class Data> |
| 548 | Sample fractional(const Data &data, Sample fractional) const { |
| 549 | Sample subSampleDelay = fractional*subSampleSteps; |
| 550 | int lowIndex = subSampleDelay; |
| 551 | if (lowIndex >= subSampleSteps) lowIndex = subSampleSteps - 1; |
| 552 | Sample subSampleFractional = subSampleDelay - lowIndex; |
| 553 | int highIndex = lowIndex + 1; |
| 554 | |
| 555 | Sample sumLow = 0, sumHigh = 0; |
| 556 | const Sample *coeffLow = coefficients.data() + lowIndex*n; |
| 557 | const Sample *coeffHigh = coefficients.data() + highIndex*n; |
| 558 | for (int i = 0; i < n; ++i) { |
| 559 | sumLow += data[i]*coeffLow[i]; |
| 560 | sumHigh += data[i]*coeffHigh[i]; |
| 561 | } |
| 562 | return sumLow + (sumHigh - sumLow)*subSampleFractional; |
| 563 | } |
| 564 | }; |
| 565 | |
| 566 | template<typename Sample> |
| 567 | using InterpolatorKaiserSinc20 = InterpolatorKaiserSincN<Sample, 20>; |
| 568 | template<typename Sample> |
| 569 | using InterpolatorKaiserSinc8 = InterpolatorKaiserSincN<Sample, 8>; |
| 570 | template<typename Sample> |
| 571 | using InterpolatorKaiserSinc4 = InterpolatorKaiserSincN<Sample, 4>; |
| 572 | |
| 573 | template<typename Sample> |
| 574 | using InterpolatorKaiserSinc20Min = InterpolatorKaiserSincN<Sample, 20, true>; |
| 575 | template<typename Sample> |
| 576 | using InterpolatorKaiserSinc8Min = InterpolatorKaiserSincN<Sample, 8, true>; |
| 577 | template<typename Sample> |
| 578 | using InterpolatorKaiserSinc4Min = InterpolatorKaiserSincN<Sample, 4, true>; |
| 579 | /// @} |
| 580 | |
| 581 | /** @brief A delay-line reader which uses an external buffer |
| 582 | |
| 583 | This is useful if you have multiple delay-lines reading from the same buffer. |
| 584 | */ |
| 585 | template<class Sample, template<typename> class Interpolator=InterpolatorLinear> |
| 586 | class Reader : public Interpolator<Sample> /* so we can get the empty-base-class optimisation */ { |
| 587 | using Super = Interpolator<Sample>; |
| 588 | public: |
| 589 | Reader () {} |
| 590 | /// Pass in a configured interpolator |
| 591 | Reader (const Interpolator<Sample> &interpolator) : Super(interpolator) {} |
| 592 | |
| 593 | template<typename Buffer> |
| 594 | Sample read(const Buffer &buffer, Sample delaySamples) const { |
| 595 | int startIndex = delaySamples; |
| 596 | Sample remainder = delaySamples - startIndex; |
| 597 | |
| 598 | // Delay buffers use negative indices, but interpolators use positive ones |
| 599 | using View = decltype(buffer - startIndex); |
| 600 | struct Flipped { |
| 601 | View view; |
| 602 | Sample operator [](int i) const { |
| 603 | return view[-i]; |
| 604 | } |
| 605 | }; |
| 606 | return Super::fractional(Flipped{buffer - startIndex}, remainder); |
| 607 | } |
| 608 | }; |
| 609 | |
| 610 | /** @brief A single-channel delay-line containing its own buffer.*/ |
| 611 | template<class Sample, template<typename> class Interpolator=InterpolatorLinear> |
| 612 | class Delay : private Reader<Sample, Interpolator> { |
| 613 | using Super = Reader<Sample, Interpolator>; |
| 614 | Buffer<Sample> buffer; |
| 615 | public: |
| 616 | static constexpr Sample latency = Super::latency; |
| 617 | |
| 618 | Delay(int capacity=0) : buffer(1 + capacity + Super::inputLength) {} |
| 619 | /// Pass in a configured interpolator |
| 620 | Delay(const Interpolator<Sample> &interp, int capacity=0) : Super(interp), buffer(1 + capacity + Super::inputLength) {} |
| 621 | |
| 622 | void reset(Sample value=Sample()) { |
| 623 | buffer.reset(value); |
| 624 | } |
| 625 | void resize(int minCapacity, Sample value=Sample()) { |
| 626 | buffer.resize(minCapacity + Super::inputLength, value); |
| 627 | } |
| 628 | |
| 629 | /** Read a sample from `delaySamples` >= 0 in the past. |
| 630 | The interpolator may add its own latency on top of this (see `Delay::latency`). The default interpolation (linear) has 0 latency. |
| 631 | */ |
| 632 | Sample read(Sample delaySamples) const { |
| 633 | return Super::read(buffer, delaySamples); |
| 634 | } |
| 635 | /// Writes a sample. Returns the same object, so that you can say `delay.write(v).read(delay)`. |
| 636 | Delay & write(Sample value) { |
| 637 | ++buffer; |
| 638 | buffer[0] = value; |
| 639 | return *this; |
| 640 | } |
| 641 | }; |
| 642 | |
| 643 | /** @brief A multi-channel delay-line with its own buffer. */ |
| 644 | template<class Sample, template<typename> class Interpolator=InterpolatorLinear> |
| 645 | class MultiDelay : private Reader<Sample, Interpolator> { |
| 646 | using Super = Reader<Sample, Interpolator>; |
| 647 | int channels; |
| 648 | MultiBuffer<Sample> multiBuffer; |
| 649 | public: |
| 650 | static constexpr Sample latency = Super::latency; |
| 651 | |
| 652 | MultiDelay(int channels=0, int capacity=0) : channels(channels), multiBuffer(channels, 1 + capacity + Super::inputLength) {} |
| 653 | |
| 654 | void reset(Sample value=Sample()) { |
| 655 | multiBuffer.reset(value); |
| 656 | } |
| 657 | void resize(int nChannels, int capacity, Sample value=Sample()) { |
| 658 | channels = nChannels; |
| 659 | multiBuffer.resize(channels, capacity + Super::inputLength, value); |
| 660 | } |
| 661 | |
| 662 | /// A single-channel delay-line view, similar to a `const Delay` |
| 663 | struct ChannelView { |
| 664 | static constexpr Sample latency = Super::latency; |
| 665 | |
| 666 | const Super &reader; |
| 667 | typename MultiBuffer<Sample>::ConstChannel channel; |
| 668 | |
| 669 | Sample read(Sample delaySamples) const { |
| 670 | return reader.read(channel, delaySamples); |
| 671 | } |
| 672 | }; |
| 673 | ChannelView operator [](int channel) const { |
| 674 | return ChannelView{*this, multiBuffer[channel]}; |
| 675 | } |
| 676 | |
| 677 | /// A multi-channel result, lazily calculating samples |
| 678 | struct DelayView { |
| 679 | Super &reader; |
| 680 | typename MultiBuffer<Sample>::ConstView view; |
| 681 | Sample delaySamples; |
| 682 | |
| 683 | // Calculate samples on-the-fly |
| 684 | Sample operator [](int c) const { |
| 685 | return reader.read(view[c], delaySamples); |
| 686 | } |
| 687 | }; |
| 688 | DelayView read(Sample delaySamples) { |
| 689 | return DelayView{*this, multiBuffer.constView(), delaySamples}; |
| 690 | } |
| 691 | /// Reads into the provided output structure |
| 692 | template<class Output> |
| 693 | void read(Sample delaySamples, Output &output) { |
| 694 | for (int c = 0; c < channels; ++c) { |
| 695 | output[c] = Super::read(multiBuffer[c], delaySamples); |
| 696 | } |
| 697 | } |
| 698 | /// Reads separate delays for each channel |
| 699 | template<class Delays, class Output> |
| 700 | void readMulti(const Delays &delays, Output &output) { |
| 701 | for (int c = 0; c < channels; ++c) { |
| 702 | output[c] = Super::read(multiBuffer[c], delays[c]); |
| 703 | } |
| 704 | } |
| 705 | template<class Data> |
| 706 | MultiDelay & write(const Data &data) { |
| 707 | ++multiBuffer; |
| 708 | for (int c = 0; c < channels; ++c) { |
| 709 | multiBuffer[c][0] = data[c]; |
| 710 | } |
| 711 | return *this; |
| 712 | } |
| 713 | }; |
| 714 | |
| 715 | /** @} */ |
| 716 | }} // signalsmith::delay:: |
| 717 | #endif // include guard |
| 718 | |