Polynominal Waveshaper

notes
The following code will describe how to excite discrete harmonics and only these
harmonics. A simple polynominal waveshaper for processing the data is included as well.
However the code don't claim to be optimized. Using a horner scheme with precalculated
coefficients should be your choice here.
Also remember to oversample the data (optimal in the order of the harmonics) to have them
alias free.
code
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
We assume the input is a sinewave (works for any input signal, but this makes everything more clear).
Then we have x = sin(a)

the first harmonic is plain simple (using trigonometric identities):

cos(2*a)= cos^2(a) - sin^2(a) = 1 - 2 sin^2(a)

using the general trigonometric identities:

sin(x + y) = sin(x)*cos(y) + sin(y)*cos(x)
cos(x + y) = cos(x)*cos(y) - sin(y)*sin(x)

together with some math, you can easily calculate: sin(3x), cos(4x), sin(5x), and so on...


Here's how the resulting waveshaper may look like:

// o = output, i = input
o = fPhase[1]*     i                                                                  * fGains[0]+
    fPhase[1]*(  2*i*i             -   1                                            ) * fGains[1]+
    fPhase[2]*(  4*i*i*i           -   3*i                                          ) * fGains[2]+
    fPhase[3]*(  8*i*i*i*i         -   8*i*i         +   1                          ) * fGains[3]-
    fPhase[4]*( 16*i*i*i*i*i       -  20*i*i*i       +   5 * i                      ) * fGains[4]+
    fPhase[5]*( 32*i*i*i*i*i*i     -  48*i*i*i*i     +  18 * i*i     -  1           ) * fGains[5]-
    fPhase[6]*( 64*i*i*i*i*i*i*i   - 112*i*i*i*i*i   +  56 * i*i*i   -  7 * i       ) * fGains[6]+
    fPhase[7]*(128*i*i*i*i*i*i*i*i - 256*i*i*i*i*i*i + 160 * i*i*i*i - 32 * i*i + 1 ) * fGains[7];

fPhase[..] is the sign array and fGains[..] is the gain factor array.

P.S.: I don't want to see a single comment about the fact that the code above is unoptimized. I know that!

Comments

Here's the more math like version:

// o = output, i = input
o = fPhase[1]* i * fGains[0]+
    fPhase[1]*(  2*i^2 - 1 ) * fGains[1]+
    fPhase[2]*(  4*i^3 - 3*i ) * fGains[2]+
    fPhase[3]*(  8*i^4 - 8*i^2 + 1 ) * fGains[3]-
    fPhase[4]*( 16*i^5 - 20*i^3 + 5*i ) * fGains[4]+
    fPhase[5]*( 32*i^6 - 48*i^4 + 18 * i^2 - 1 ) * fGains[5]-
    fPhase[6]*( 64*i^7 - 112*i^5 + 56 * i^3 - 7 * i ) * fGains[6]+
    fPhase[7]*(128*i^8 - 256*i^6 + 160 * i^4 - 32 * i^2 + 1 ) * fGains[7];

Actually, this *doesn't* work in the way described on any input, only on single sinusoids of
amplitude 1. It's nonlinear - (a+b)^n is not the same thing as a^n+b^n, nor are (a*b)^n and a*(b^n).
Even just a sum of two sinusoids or a single sinusoid of a different amplitude breaks the
chosen-harmonics-only thing.

Do'oh. You're right. Once more I got fooled by the way of my measurement. That explains a
lot of things... Thanks for the clarification!

Btw. the coeffitients follow the chebyshev polynomials. Just in case you wonder about the
logic behind. Maybe we can call it chebyshev waveshaper from now on...

I played with this idea for a while yesterday to no avail before discovering this post tonight.
I thought I could excite any harmonic I wanted using select Chebyshev Polynomials. But you are
totally right - it doesn't work that way. Any complex waveform that can be broken down into a
Fourier series is a linear sum of terms. Squaring or cubing the waveform, and therefor this sum,
leads to multiple cross terms which introduce additional frequencies. It does only work with
normalized single sinusoids . . .which is too bad.

Right now, the only way I can see to do this sort of thing where you excite select harmonics at will
is to run an FFT and then work from there in the frequency domain.

But my question is, if we are looking to simulate tube saturation, is the Chebyshev method good
enough. What, after all, do tubes do? Does a tube amp actually add discrete harmonics or is it
introducing all of those cross term frequencies as well?

According to another post, the tube is simply a non linear function, for example a tan(x).
Actually by saturating any signal you will get harmonics (any but a pure square which cannot
be more saturated of couse...). As tan(x) is not linear, you should get harmonics… that's all.
Now if you want to pass only the high frequencies,just split the signal into 2 frequencies using
a lowpass vs highpass = signal - lowpass and process the frequencies you want to.

I would like to add that this method could result in an approximated harmonic exciter using an
array of filters sufficiently narrow and steep to approximately single out individual frequencies
composing the original signal.

As such, it would be processing intensive because the polynomial would need to be calculated on
each band.

What you have presented is not completely bad.  You only need to take into consideration that you're
getting frequency terms that are not necessarily harmonics.  Steve Harris has a LADSPA plugin that
uses the chebychev polynomial as a waveshaper...and he calls it a harmonic exciter.

To user physicstrick:  Tube emulation is much more than waveshaping.  Bias conditions change with
signal dynamics, and you essentially get signal-power modulated duty cycle.  I have found some good
articles about this and also there is a commercial product that claims to solve the discretized
system of ODE's in real time.  This model eats CPU like you would not imagine.

My simple "trick" is to include the nonlinear function in a 1rst order filter calculation and
also to modulate the filter time constants with the filter state variable amplitude.  This is not
quite right, but it produces an emulation that is more pleasing than plain waveshaping.

If I'm correct, you're describing the same technique achievable by use of Chebyschev polynomials,
i.e. generating any integral harmonic of the original signal. I've experimented with these,
synthesizing only the second, third, fourth, etc. harmonics, but could never get a realistic
sound, probably because overdiven tubes/transistors don't work this way and there's no intermodulation
distortion, only the pure harmonics.