// True-peak L-upsample windowed-sinc, long double. Not ITU BS.1770.
// g++ -O2 -std=c++17 -o /tmp/tp_ld /tmp/tp_ld.cpp
#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
static constexpr long double PI = 3.141592653589793238462643383279502884L;
static long double sinc(long double x) {
  if (fabsl(x) < 1e-18L) return 1.0L;
  return sinl(PI * x) / (PI * x);
}
static long double blackman(int n, int N) {
  if (N <= 1) return 1.0L;
  long double i = (long double)n / (long double)(N - 1);
  return 0.42L - 0.5L * cosl(2 * PI * i) + 0.08L * cosl(4 * PI * i);
}
static double peak_db(const std::vector<long double> &x) {
  long double m = 0;
  for (auto v : x) m = std::max(m, fabsl(v));
  return (double)(20.0L * log10l(std::max(m, 1e-18L)));
}
static std::vector<long double> upsample(const std::vector<long double> &x, int L, int taps) {
  int N = L * taps;
  std::vector<long double> h(N);
  long double sum = 0;
  int mid = (N - 1) / 2;
  for (int n = 0; n < N; n++) {
    long double xx = (n - mid) / (long double)L;
    h[n] = sinc(xx) * blackman(n, N);
    sum += h[n];
  }
  for (auto &v : h) v *= (long double)L / sum;
  std::vector<long double> y(x.size() * (size_t)L, 0);
  int htaps = taps;
  for (size_t i = 0; i < x.size(); i++) {
    for (int p = 0; p < L; p++) {
      long double acc = 0;
      for (int k = 0; k < htaps; k++) {
        int j = (int)i - k;
        if (j >= 0) acc += h[k * L + p] * x[j];
      }
      y[i * L + p] = acc;
    }
  }
  return y;
}
int main() {
  const int n = 256;
  std::vector<long double> sine(n), clip12(n), clip48(n);
  for (int i = 0; i < n; i++) {
    sine[i] = 0.5L * sinl(2 * PI * 0.25L * i + PI / 4);
    long double s = 0.92L * sinl(2 * PI * 0.25L * i);
    clip12[i] = tanhl(12.0L * s);
    clip48[i] = tanhl(48.0L * s);
  }
  printf("input_sine_db %.4f\n", peak_db(sine));
  int Ls[] = {4, 8, 12, 32, 48};
  int taps = 16;
  for (int L : Ls) {
    auto y4 = upsample(sine, L, taps);
    auto y12 = upsample(clip12, L, taps);
    auto y48 = upsample(clip48, L, taps);
    printf("L=%d taps=%d sine_tp %.4f clip12_tp %.4f clip48_tp %.4f\n",
           L, taps, peak_db(y4), peak_db(y12), peak_db(y48));
  }
  return 0;
}