True Peak detection II

Question

As a followup to my True Peak detection question, I'm trying to implement a detection method by following this documentation using the Catmull-Rom interpolation method.

What I've done so far can be seen in my Desmos worksheet (includes sample of y -data which creates the issue in question).

I've (probably) got stationary points solving to work correctly but, calculation used with the cases not having critical points and are solved through finding inflection point using equation $x = \frac{-b}{3a}$ are not giving proper results. I suppose the result could be at least improved by using suitable conditional clauses. As seen in the Desmos worksheet, resulting x values are used as an input for the interpolation function to get the peak level solved (yi = interpolation(y, x), where y is up-sampled data, interpolated using half-polyphase FIR low-pass filter).

Any suggestions to get this issue solved?

Edit: Here are some C++ code at Compiler Explorer.

catmullrom4() is the function name related to this issue:

float catmullrom4(float* y){   
   // Solve tp through stationary and inflection points
   float a, b, c, d, x, yi;
a = -0.5fy[0] + 1.5fy[1] - 1.5fy[2] + 0.5fy[3];
   b = y[0] - 2.5fy[1] + 2.fy[2] - 0.5fy[3];
   c = -0.5fy[0] + 0.5f*y[2];
   d = y[1];
float s1 = (-b - std::sqrt((bb) - (3.0fac))) / (3.0fa);
  float s2 = (-b + std::sqrt((bb) - (3.0fac))) / (3.0fa);
  float s3 = -b/(3.0f*a);
float abs1 = std::abs(s1);
  float abs2 = std::abs(s2);
  float abs3 = std::abs(s3);
if(abs2 > abs1 ){
    x = s1;
  }
  else if (abs1 > abs2){
      x = s2;
  }
  else{
     x = s3;
    if (abs3 > 3.0f){ 
          x = 1.0f;
          }
    else if (abs3 < 2.0f){
        x = 0.0f;
        }
    }
  return (x * (x * (a * x + b) + c) + d);
}

EDIT: Updated the Desmos and Compiler Explorer links

Answer 1

I think you are suffering needlessly here. This is a well known problem with an equally well known solution. Below is the code for up sampling by 4 with a 32 tap FIR filter. The error for the worst case sine wave is less than 0.1dB and for the outlandish double impulse it's 0.3dB.

This can be implemented efficiently as a polyphase filter with 8 taps each. The 0th phase doesn't need to be calculated since it's a pure delay. So you need a total of 24 multiplies per sample (no division or transcendental functions). You can further optimize using symmetry properties: the 1rst and 3rd phases are time reversals of each other and the 2nd phase is the time reversal of itself. So there are only 12 unique coefficients.

If that's not good enough, please explain why.

%% find true peak to upsampling
%% test with a worst case sine wave: true peak error is 1.414
nx = 256;
x1 = sin(2pi(0:nx-1)'/4+pi/4);
% x = 0x;
% x(100:101) = 1;
y1=truePeak(x1,4,32);
fprintf('Sine Wave true peak error = %6.2fdB\n',20log10(max(y1)));
%% test with a double impulse
x2 = zeros(nx,1);
gain = 1/(2sinc(0.5));
x2(nx/2) = gain;
x2(nx/2+1) = gain;
y2=truePeak(x2,4,32);
fprintf('Double impulse true peak error = %6.2fdB\n',20log10(max(y2)));
%% true peak finder (not optimized)
function ymax = truePeak(x,nUp,nFir)
% desugn the lowpass filter as a windowed sinc
t = (-nFir/2:nFir/2)'/nUp;
h = sinc(t).*hamming(nFir+1);
nx = length(x);
y = zeros(nx,nUp);
for i=1:nUp
  y(:,i) = conv(x,h(i:nUp:end),'same');
end
ymax = max(abs(y),[],2);
end

Answer 2

If the goal is to identify inter-sample overs such that you can normalize a signal and still not have voltages out of a D/A converter exceeding some maximum, you are essentially modelling the D/A conversion process, and without that knowledge, I think that there can be no absolute limit.

The D/A could be filtered by a close approximation to a sinc - or it could not.

Answer 3

OK, followed @robertbristow-johnson's suggestion in comments and implemented TP detection based on quadratic functions.

Here's the C++ code at Compiler Explorer site with 4x up-sampled data from test case in Hilmar̈́s answer. Given C++ implementation needs some modifications to work with different up-sampling ratios.

Here's also Desmos sheet demonstrating the main functionality of the implementation.