Monotonic, symmetrical soft-clipping polynomial

Question

Okay, I intend to answer this question myself, but I would like other folks to pipe in on it. The issues or requirements are this:

1. Input to soft-clipper is $x(t)$, output is $y(t)$.
2. Output $y(t)$ limited in magnitude to $1$: $|y(t)| \le 1$
3. Memoryless. $y(t)$ only a function of the present $x(t)$, no past values. So we may as well call the input $x$ and the output $y$. And we will call the memoryless soft-clipping function: $ y \triangleq f(x) $.
4. $f(-1)=-1$ and $f(1) = 1$.
5. $f(-x) = -f(x)$, an odd-symmetry (about $x=0$) function (so polarity of the input is irrelevant).
6. Monotonically increasing function: $f'(x) \ge 0 \quad \forall x$
7. For $|x| \le 1$, then $f(x)$ is a finite-order polynomial (so the frequency components generated are finite in frequency) and an odd-order polynomial (so all even-order terms are zero because of the odd-symmetry). Let's call that odd order $2N+1$ where $N\in \mathbb{Z} \ge 0$. Being a finite-order polynomial means a limit to the width of spurious frequency components generated which tells us what oversampling ratio suffices.
8. $f(x)$ to have as many derivatives equal to zero as possible where $x=\pm 1$. So at $\pm 1$, this splices to a constant function ($f(x) = \operatorname{sgn}(x)$ for $|x| \ge 1$) with as many derivatives being continuous as possible.

So what is $f(x)$?

$$ f(x) = \begin{cases} -1 & x \le -1 \\ \sum\limits_{n=0}^{N} a_n \, x^{2n+1} \quad & -1 \le x \le +1 \\ +1 & +1 \le x \\ \end{cases}$$

$$ $$

$$\begin{align} y(t) & \triangleq f\big(x(t)\big) \\ & = \sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n+1} \\ & = x(t)\sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n} \quad |x(t)| \le 1 \end{align}$$

$f(x)$ is continuous everywhere and as many derivatives as possible are continuous everywhere and the only discontinuity in any derivative is at $x = \pm 1$.

What are the odd-order polynomial coefficients $a_n$?

Answer 1

Here's a recursive formula for the soft-clipping polynomials $f_N(x)$ of degree $2N+1$:

$$\begin{align}f_0(x) &= x\\ f_N(x) &= f_{N-1}(x) + \frac{(2N)!}{4^{N}(N!)^2}\,\left(1-x^2\right)^N x\\ &= \sum_{n=0}^{N}\frac{(2n)!}{4^{n}(n!)^2}\,\left(1-x^2\right)^n x\end{align}$$

with splice points $f_N(\pm1) = \pm1$ at each of which the first $N$ derivatives equal zero. The clunky rational number coefficients are equal to OEIS A001790 divided by OEIS A046161 from the On-Line Encyclopedia of Integer Sequences (OEIS). Maybe it's cheating but whenever I see the beginnings of a sequence of integer numbers in my math formulas, I go there to search and they virtually always have a formula for the numbers. Rational numbers are usually split into a numerator and a denominator OEIS entry, not sure if with common factors.

If you write $\sin(x)$ in place of $x$ in the function body, you can use these as band-limited square wave approximations free of overshoot and ripple. I think that as $N\to\infty$, the coefficients of the expanded polynomial approach the coefficients of the Fourier series of a square wave. (At least the coefficient of $x$ approaches $\frac{4}{\pi}$, which is the coefficient of the fundamental frequency in the Fourier series.)

If the slope at $x=0$ is normalized to $1$ by $g(x) = f\left(\frac{x}{f'(0)}\right)$, where $f'(x)$ is the derivative of $f(x),$ then $g_\infty(x)$ seems to approach the integral of a Gaussian function:

Figure 1. $g_{75}(x)$ (blue) and $\operatorname{ERF}\left(\frac{\sqrt{\pi}x}{2}\right)$ (red) with identical slope at $x=0$.

This is confirmed by that we can stretch horizontally the integrand $(1-u^2)^N$ (see Robert's answer) by a substitution $u \to \frac{u}{\sqrt{N}}$ so that in the limit $N\to\infty$ it becomes Gaussian:

$$\lim_{N\to\infty}\left(1 - \left(\frac{u}{\sqrt{N}}\right)^2\right)^N = e^{-u^2}.$$

Answer 2

Okay, so to satisfy requirements 5, 6, 7, and 8, the soft-clipping function will have the form:

$$f(x) = \begin{cases} -1 & x \le -1 \\ K \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du \quad & -1\le x\le 1 \\ +1 & 1 \le x \\ \end{cases}$$

where $K$ is a constant judiciously chosen to satisfy requirement 4.

$$ K = \left( \int\limits_{0}^{1} \big(1 - u^2 \big)^N \ du \right)^{-1} $$

The 1st derivative of $f(x)$ is

$$ f'(x) = K \big(1 - x^2 \big)^N $$

The 2nd derivative of $f(x)$ is

$$ f''(x) = K N \big(1 - x^2 \big)^{N-1} (-2x)$$

The 3rd derivative of $f(x)$ is

$$\begin{align} f'''(x) &= K N(N-1) \big(1 - x^2 \big)^{N-2} (-2x)^2 \ + \ K N \big(1 - x^2 \big)^{N-1} (-2) \\ &= K \big(1 - x^2 \big)^{N-2} \bigg( N(N-1)(-2x)^2 - 2N \big(1 - x^2 \big) \bigg) \\ \end{align}$$

and, for $n \ge 1$, the $n$th derivative is

$$ f^{(n)}(x) = K \big(1 - x^2 \big)^{N-n+1} \, g_n(x) $$

where $g_n(x)$ is some ($n$-1)th order polynomial function of $x$ and is finite in value. This can be proven inductively by considering the ($n$+1)th derivative:

$$\begin{align} \\ f^{(n+1)}(x) &= \tfrac{d}{dx} \Big( f^{(n)}(x) \Big) \\ &= \tfrac{d}{dx} \Big(K \big(1 - x^2 \big)^{N-n+1} \, g_n(x)\Big) \\ &= K(N-n+1)\big(1 - x^2 \big)^{N-n} (-2x) g_n(x) \, + \, K \big(1-x^2 \big)^{N-n+1}g'_n(x) \\ &= K\big(1 - x^2 \big)^{N-n} \Big( (N-n+1)(-2x) g_n(x) \, + \, (1-x^2) g'_n(x) \Big) \\ &= K\big(1 - x^2 \big)^{N-n} \, g_{n+1}(x) \\ \\ \end{align}$$

where $ g_{n+1}(x) = (N-n+1)(-2x) g_n(x) + (1-x^2) g'_n(x) $.

Because of differentiation, the order of polynomial $g'_n(x)$ is one less than the order of $g_n(x)$, but polynomial $(1-x^2)g'_n(x)$ is one order greater than $g_n(x)$ and so also is $(-2x) g_n(x)$.

When $x = \pm 1$, then the first $N$ derivatives are zero, $$ f^{(n)}(x) \Bigg|_{x=\pm 1} = 0 \qquad \text{for } 1 \le n \le N $$ making this polynomial maximally flat at $x = \pm 1$.

The integrand is a binomial and can be expressed as a power series using binomial expansion:

$$\begin{align} \big(1 - u^2 \big)^N & = \sum\limits_{n=0}^{N} \binom{N}{n} \big(-u^2\big)^n (1)^{N-n} \\ & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \big(-u^2\big)^n (1)^{N-n} \\ & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \\ \end{align}$$

So the integral can be expressed as an integral of a power series:

$$\begin{align} \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du & = \int\limits_{0}^{x} \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \ du \\ & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n \int\limits_{0}^{x} u^{2n} \ du \\ & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} u^{2n+1} \Bigg|_0^x \\ & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} x^{2n+1} \\ & = x \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \big(x^2 \big)^n \\ \end{align}$$

When $x = \pm 1$, we get

$$\begin{align} \int\limits_{0}^{\pm 1} \big(1 - u^2 \big)^N \ du &= \pm \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \\ &= \pm \frac{4^N}{2N+1}\binom{2N}{N}^{-1}\\ \end{align}$$

(there is a little improperness in the integral limit.) The bottom equality is because of this Appendix proof. So the scaler $K$ must be

$$\begin{align} K &= \left( \frac{4^N}{2N+1}\binom{2N}{N}^{-1} \right)^{-1} \\ \\ &= \frac{2N+1}{4^N}\binom{2N}{N} \\ \\ &= \frac{2N+1}{4^N} \frac{(2N)!}{(N!)^2} \\ \\ &= \frac{(2N+1)!}{4^N (N!)^2} \\ \end{align}$$

This makes the entire soft-clipping function to be:

$$ f(x) = \begin{cases} -1 & x \le -1 \\ \\ x \cdot \frac{(2N+1)!}{4^N \ N!} \sum\limits_{n=0}^{N} \frac{(-1)^n}{(2n+1)\,n!\,(N-n)!} \big(x^2 \big)^n \quad & -1 \le x \le +1 \\ \\ +1 & +1 \le x \\ \end{cases} $$

and it appears that the odd-power (that is the power = $2n+1$) polynomial coefficients are

$$ a_n = (-1)^n \frac{(2N+1)!}{4^N \, N! \,(2n+1) \, n! \, (N-n)!} $$

The polynomials (without splicing to the $\operatorname{sgn}(x) = \pm 1$ saturated components) look like

I think the order $2N+1$ starts at 1 and goes to 9 (or $0 \le N \le 4$)

With the saturation attached, the curves look like

The soft-clipping function is continuous everywhere and all derivatives, up to the $N$th derivative are continuous everywhere and the ($N$+1)th derivative and higher are continuous everywhere except at the splices where $x = \pm 1$.

Here are the same set of curves but with the scaling adjusted so that the slope around $x=0$ (or "gain") remains at 1 (or "0 dB gain").

The point of discontinuity (where the polynomial is spliced to a constant $\pm 1$) is at $x= \pm K$. I really don't think one needs to get over a $7$th-order ($N=3$) softclipper. Oversampling by a factor $4$x suffices to deal with this $7$th-order polynomial and prevent any aliasing to the original baseband before upsampling.

Answer 3

There is a similar group of functions called smoothstep:

https://en.wikipedia.org/wiki/Smoothstep

Apart from the different range, they seem to be identical to they above, but could be easier to calculate because they have less nonzero coefficients.

Answer 4

Cool stuff! Not really an answer but a few comments, if I may

A few years ago, I did an exercise trying to create a "fade in" polynomial window that pase zero derivatives up to order $n$ at both edges. That feels like a similar exercise to yours with some substitution $$f(x) = c \cdot g(a \cdot x + b) + d $$ or specifically $f(x) = 2 \cdot g(2 \cdot x - 1) -1$.

The derivative of such a polynomial needs to have shape of $$\frac{\partial }{\partial x}g(x)=x^n \cdot (x-1)^n$$ You can start with a polynomial like this, integrate once, and normalize to $g(1)=1$. Matlab code below.

I never bothered to derive a closed form solution, though because it's easy enough to calculated them the way I described above and it was less useful than I had hoped

1. The higher you drive the order, the more it resembles a step function or hard clipper respectively. Helps a bit with high frequencies but is bad for low frequencies.
2. It does not make a great soft clipper for audio. I found it much more useful to do soft clipping by stitching together a linear 1:1 section with a 2nd-4th order polynomial that flattens it out. That poly needs something like $g(x_{lin}) = x_{lin}; g^{'}(x_{lin}) = - 1; g(x_{clip}) = 1; g^{'}(x_{clip}) = 0$

function p = polwin3(n)
% creates a polynomial half-window which which has 0 derivatives up
% to order n. p(0) = 0; p(1) = 1;

% this version actually constructs the polynomial from the derivative 
% polynomial which x^n * (x-1)^n

% calculate top  half of derivative poly (x-1)^n
% lower half is all zero anyway
b = binom(n);
% alternate the sign since it's (x-1)
b(2:2:length(b)) = -1 * b(2:2:length(b)) ;
% integrate
b = b ./ ((2*n+1):-1:(n+1));
% build plynomial, add the trailing zeros
p = [b zeros(1,n+1)];
% normalize to p(1) = 1;
p = p ./ polyval(p,1);
% go to integer
p = round(p);
end

function b = binom(n)
% calculate binomial coefficients for order n
b = [1 1];
for i = 2:n
  b = [1 b(1:end-1)+b(2:end) 1];
end