# Monotonic, symmetrical soft-clipping polynomial

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$?

• Hi! Great Question (And answer), but would you mind adding a reference for the statement between parenthesis in point 7? it is not immediately obvious that an odd order polynomial clipper would generate finite harmonics and all odd order as well. Maybe there is a simple proof for this but would like to see it :) – bone Apr 24 '17 at 9:58
• hi @bone , so i split the parenthed remarks into two parenthed remarks to associate them more directly with their associated statements. (a sinusoid that is processed by an $N$-th order polynomial will generate, at highest, a sinusoid of $N$ times the frequency.) – robert bristow-johnson Apr 24 '17 at 10:03
• Why a polynomial? Erf and logistic are common soft clippers – Stanley Pawlukiewicz Apr 8 '18 at 16:25
• well i guess, @StanleyPawlukiewicz, the main reason is, at least when the signal output is not saturated, then "the frequency components generated [by the non-linearity] are finite in frequency" if you know it's limited to 7th-order, you will have no components higher than 7× Nyquist coming out and you can manage that with 4× oversampling (the top 3 generated overtones fold over, but do not contaminate the original baseband). – robert bristow-johnson Apr 9 '18 at 4:27

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.

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

$$\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}$$

(there is a little improperness in the integral limit.) So the scaler $K$ must be

$$K = \left( \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \right)^{-1}$$

This makes the entire soft-clipping function to be:

$$f(x) = \begin{cases} -1 & x \le -1 \\ \\ \left( \sum\limits_{i=0}^{N} \frac{N!}{i!(N-i)!} \frac{(-1)^i}{2i+1} \right)^{-1} x \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \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 = \left( \sum\limits_{i=0}^{N} \frac{N!}{i!(N-i)!} \frac{(-1)^i}{2i+1} \right)^{-1} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1}$$

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.

• now if someone can add an answer that turns these polynomials into a sum of Tchebyshev polynomials, that would be maybe useful. – robert bristow-johnson Dec 11 '16 at 7:05
• it got spurred on by this SE post and then i had an email exchange with Viktor Toth about his method to avoid clipping: $$y(t) = x_1(t) + x_2(t) - x_1(t)x_2(t)$$ – robert bristow-johnson Dec 12 '16 at 2:00
• that smooth complementary splicing is a very similar problem to the complementary highpass/lowpass filter banks problem. like in this little thing. – robert bristow-johnson Dec 12 '16 at 2:26
• Turning $f(x)$ into Tchebyshev polynomials is do-able, but moderately hideous. From what i worked out we have $$f(x) = \sum_{n=0}^N b_{k} T_{2k+1}$$ where $$b_k = \sum_{n=0}^{N-k}\left(\begin{array}2(n+k)\\n\end{array}\right)\frac{a_{2(n+k)+1)}}{2^n}$$ where $a_n$ are the coefficients you already calculated. – Peter Cudmore Dec 13 '16 at 12:35
• @robertbristow-johnson No worries, it make take some time since i no longer have my notes. I'll post the derivation as an answer. Also, whoops. – Peter Cudmore Aug 1 '17 at 23:48

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}.$$

• so does this satisfy $$f'_N(x) = K_N \big(1 - x^2 \big)^N$$ with $$K_N = \frac{1}{ \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} }$$ ? i think you mean that it does, Olli.  BTW, i think you can see how replacing $x$ with $\cos(\omega)$ yields that same "square" approximation (without ripples) in the frequency domain and it's time-limited, an FIR. i think the wavelet and filterbank people make use of these for to get both perfect reconstruction and to separate the frequency bands better and better. – robert bristow-johnson Apr 17 '17 at 23:30
• Yes, these are the same as your polynomials. You need to use sines for the square wave, as using cosines does a Hilbert transform. Yes I can see that it would work as a filter. – Olli Niemitalo Apr 18 '17 at 5:51
• WTF?? @Olli. that makes no fucking sense. as a memoryless, time-invariant soft-clipping polynomial mapping, whether you drive your input with $\sin(\omega t)$ or with $\sin(\omega (t - \tfrac{-\pi}{2\omega}))$ should make no difference except a delay of $\frac{-\pi}{2\omega}$ in the output. – robert bristow-johnson Apr 18 '17 at 18:51
• You're absolutely right. I got confused because I did $x \to \sin(x)$, and expanded the formula to have harmonic sines and then did $\sin(kx) \to \cos(kx)$. But that's not the same as going $x \to \cos(x)$ directly. – Olli Niemitalo Apr 18 '17 at 19:11

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.

• hey, Stefan. the wikipedia page seems to have problems, but the microsoft reference says that they're Hermite. if smoothstep is Hermite polynomials, i suspect they are the same polynomials as i have derived. – robert bristow-johnson Dec 18 '16 at 19:45
• i have asked a question at the Math SE to try to nail down the relationship between this and Smoothstep. – robert bristow-johnson Apr 24 '17 at 6:26
• just to let you know Stefan, that i had not previously known of this Smoothstep function definition. we have picked up this relation and have proven that this soft-clipping polynomial with odd-symmetry is, indeed, the Smootstep function with a translation and scaling in both the $x$-axis and $y$-axis. we were able to prove the relationship explicitly. and i have started modifying Wikipedia a little to reflect this relationship. the number of nonzero coefficients should be the same.thank you for bringing this to my attention and connecting this. – robert bristow-johnson May 3 '17 at 4:33

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