You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

83 lines
1.6KB

  1. #include <dsp/minblep.hpp>
  2. #include <dsp/fft.hpp>
  3. #include <dsp/window.hpp>
  4. namespace rack {
  5. namespace dsp {
  6. void minBlepImpulse(int z, int o, float* output) {
  7. // Symmetric sinc array with `z` zero-crossings on each side
  8. int n = 2 * z * o;
  9. float* x = new float[n];
  10. for (int i = 0; i < n; i++) {
  11. float p = math::rescale((float) i, 0.f, (float)(n - 1), (float) - z, (float) z);
  12. x[i] = sinc(p);
  13. }
  14. // Apply window
  15. blackmanHarrisWindow(x, n);
  16. // Real cepstrum
  17. float* fx = new float[2 * n];
  18. RealFFT rfft(n);
  19. rfft.rfft(x, fx);
  20. // fx = log(abs(fx))
  21. fx[0] = std::log(std::fabs(fx[0]));
  22. for (int i = 1; i < n; i++) {
  23. fx[2 * i] = std::log(std::hypot(fx[2 * i], fx[2 * i + 1]));
  24. fx[2 * i + 1] = 0.f;
  25. }
  26. fx[1] = std::log(std::fabs(fx[1]));
  27. // Clamp values in case we have -inf
  28. for (int i = 0; i < 2 * n; i++) {
  29. fx[i] = std::fmax(-30.f, fx[i]);
  30. }
  31. rfft.irfft(fx, x);
  32. rfft.scale(x);
  33. // Minimum-phase reconstruction
  34. for (int i = 1; i < n / 2; i++) {
  35. x[i] *= 2.f;
  36. }
  37. for (int i = (n + 1) / 2; i < n; i++) {
  38. x[i] = 0.f;
  39. }
  40. rfft.rfft(x, fx);
  41. // fx = exp(fx)
  42. fx[0] = std::exp(fx[0]);
  43. for (int i = 1; i < n; i++) {
  44. float re = std::exp(fx[2 * i]);
  45. float im = fx[2 * i + 1];
  46. fx[2 * i] = re * std::cos(im);
  47. fx[2 * i + 1] = re * std::sin(im);
  48. }
  49. fx[1] = std::exp(fx[1]);
  50. rfft.irfft(fx, x);
  51. rfft.scale(x);
  52. // Integrate
  53. float total = 0.f;
  54. for (int i = 0; i < n; i++) {
  55. total += x[i];
  56. x[i] = total;
  57. }
  58. // Normalize
  59. float norm = 1.f / x[n - 1];
  60. for (int i = 0; i < n; i++) {
  61. x[i] *= norm;
  62. }
  63. std::memcpy(output, x, n * sizeof(float));
  64. // Cleanup
  65. delete[] x;
  66. delete[] fx;
  67. }
  68. } // namespace dsp
  69. } // namespace rack