The method is based on previous work identifying simple closed form equations for a finite sum of harmonics with constant or exponential decaying coefficients (e.g. a band limited impulse train "BLIT"), and using that to closely approximate inverse linear and inverse quadratically decaying coefficients (e.g. triangle, square, and saw waves).
The motivation for this lies in the weaknesses of alternative methods for generating such waves. Methods relying on integration struggle to stay centered under frequency modulation as there is interesting transient behavior. Wavetables can require a lot of memory and have a trade-off between interpolation kernel performance and alias-limiting. Up-sampling and filtering methods have a similar tradeoff and both can sound muffled as they tend to roll-off early at the high end of the spectrum. And simple additive synthesis even with trig identity optimizations becomes inefficient compared to this method above 16 harmonics, which is inadequate fidelity for most applications. This method aims to address all these drawbacks by providing a closed form computation per sample which has constant-time performance, no transient instability, and absolutely zero aliasing for constant frequency input.
The previously known closed form equation for BLIT functions, with f (frequency in Hertz), t (time in seconds), and n (number of harmonics, generally (sample rate)/(2f)) is:
$$ \text{Let} \quad \theta = (2\pi f \cdot t + \phi) $$ $$ \text{blit}_{n}(\theta) = \sum_{k=0}^{n} \frac{1}{n} \cos(k\theta) = \frac{\sin((2n+1)\frac{\theta}{2})}{2n \cdot \sin(\frac{\theta}{2})} $$Here is a graph of that function at n=15
Now that we have established this function, let's look at the Fourier series for a band limited sawtooth wave:
$$ \text{saw}_n(\theta) = \frac{1}{\pi} \sum_{k=1}^{n} \frac{1}{k} \sin(k\theta) $$Here is a plot of 1/k (in red), five A exp(-Bk) functions (in blue), and their sum (in green), up to k=100 on a logarithmic y axis:
That tells us we can approximate the sawtooth function very accurately like so:
$$ \text{saw}_n(\theta) = \frac{1}{\pi} \sum_{k=1}^{n} \frac{1}{k} \sin(k\theta) \approx A_1 \sum_{k=1}^{n} e^{-B_1 k} \sin(k\theta) + A_2 \sum_{k=1}^{n} e^{-B_2 k} \sin(k\theta) + ...$$Finding the closed form calculation of that new sum...
$$ A \sum_{k=1}^{n} e^{-Bk} \sin(k\theta) $$ $$ = A \cdot \text{Imag} \left\{ \sum_{k=1}^{n} e^{-Bk} e^{ki\theta} \right\} $$ $$ = A \cdot \text{Imag} \left\{ \sum_{k=1}^{n} e^{(-B+i\theta)^{k}} \right\} $$ $$ = A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - e^{(n+1)(-B+i\theta)} } { 1 - e^{-B+i\theta} } \right \} $$ $$ = A \cdot \text{Imag} \left \{ \frac{ (e^{-B+i\theta} - e^{(n+1)(-B+i\theta)})(1 - e^{(-B-i\theta)}) } { (1 - e^{-B+i\theta})(1 - e^{-B-i\theta}) } \right \} $$ $$ = A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - e^{(n+1)(-B+i\theta)} - e^{-B+i\theta} e^{-B-i\theta}) + e^{(n+1)(-B+i\theta)} e^{(-B-i\theta)}) } { 1 - e^{-B+i\theta} - e^{-B-i\theta} + e^{-2B} } \right \} $$ $$ = A \cdot \text{Imag} \left \{ \frac{ e^{-B} e^{i\theta} - e^{-(n+1)B} e^{(n+1)i\theta} - e^{-2B} + e^{-(n+2)B} e^{((n)i\theta)} } { 1 - e^{-B}(e^{i\theta} - e^{-i\theta}) + e^{-2B} } \right \} $$ $$ = A \cdot \frac{e^{-B} \sin(\theta) - e^{-(n+1)B} \sin(n\theta+\theta) + e^{-(n+2)B} \sin(n\theta) } {1-2e^{-B} cos(\theta)+e^{-2B}} $$Putting that back into our sawtooth approximation, we get this closed form sawtooth equation
$$ \text{saw}_n(\theta) = \frac{1}{\pi} \sum_{k=1}^{5} A_k \cdot \frac{e^{-B_k} \sin(\theta) - e^{-(n+1)B_k} \sin(n\theta+\theta) + e^{-(n+2)B_k} \sin(n\theta) } {1-2e^{-B_k} cos(\theta)+e^{-2B_k}}$$Computationally, it requires 4 trig functions per sample, and the exponentials only need to be recomputed when the number of harmonics changes. Additionally, the trig functions can be computed algebraically from the previous sample using the sin-cos angle-add equations.
Here is a plot of that sawtooth wave approximation using the same A and B coefficients from above, with n=35
For comparison, here is the exact version at n=35
Moving onward, let's look at the Fourier series for a band limited square wave, noting that it has the same harmonic rolloff so we can use the same coefficients as the sawtooth wave, but we need a new forumla for an odd-only cosine sum
$$ \text{square}_n(\theta) = \frac{2}{\pi} \sum_{k=1}^{\lceil n/2 \rceil} \frac{1}{2k-1} \sin((2k-1)\theta) $$The function we're looking for
$$ A \sum_{k=1}^{\lceil n/2 \rceil} e^{-B(2k-1)} \sin((2k-1)\theta) $$ $$ A \cdot \text{Imag} \left\{ \sum_{k=1}^{\lceil n/2 \rceil} e^{-B(2k-1)} e^{(2k-1)i\theta} \right\} $$ $$ A \cdot \text{Imag} \left\{ \sum_{k=1}^{\lceil n/2 \rceil} e^{(-B+i\theta)^{(2k-1)}} \right\} $$ $$ \text{Let} \quad m=\{n \quad \text{if n is odd or} \quad n-1 \quad \text{if n is even}\} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - e^{(m+2)(-B+i\theta)}} { 1 - e^{2(-B+i\theta)} } \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ (e^{-B+i\theta} - e^{(m+2)(-B+i\theta)}) \cdot (1 - e^{2(-B-i\theta)}) } { (1 - e^{2(-B+i\theta)})(1 - e^{2(-B-i\theta)}) } \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - e^{(m+2)(-B+i\theta)} - e^{-B+i\theta} e^{2(-B-i\theta)} + e^{(m+2)(-B+i\theta)} e^{2(-B-i\theta)} } { 1 - e^{2(-B+i\theta)} - e^{2(-B-i\theta)} + e^{-4B} } \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B} e^{i\theta} - e^{-(m+2)B} e^{(m+2) i\theta} - e^{-3B} e^{-i\theta} + e^{-(m+4)B} e^{(m)i\theta} } { 1 - e^{-2B}(e^{2i\theta} + e^{-2i\theta}) + e^{-4B} } \right \} $$ $$ A \cdot \frac{ e^{-B} \sin(\theta) - e^{-(m+2)B} \sin((m+2)\theta) + e^{-3B} \sin(\theta) + e^{-(m+4)B} \sin((m)\theta) } { 1 - 2e^{-2B}cos(2\theta) + e^{-4B} } $$Putting that back into our square approximation, we get (where m is odd)
$$ \text{square}_m(\theta) = \frac{2}{\pi} \sum_{k=1}^{5} A_k \cdot \frac{ (e^{-3B_k}+ e^{-B_k}) \sin(\theta) - e^{-(m+2)B_k} \sin((m+2)\theta) + e^{-(m+4)B_k} \sin(m\theta)} {1-2e^{-2B_k} cos(2\theta) + e^{-4B_k}} $$And that looks like this at n=25
For comparison, here is the exact square wave at n=25
And finally, moving on the the triangle wave:
$$ \text{triangle}_n(\theta) = \frac{4}{\pi^2} \sum_{k=1}^{\lceil n/2 \rceil} \frac{(-1)^{k-1}}{(2k-1)^2} \sin((2k-1)\theta) $$Since we're now dealing with k^-2 frequency rolloff instead of k^-1, we need a new set of exponential coefficients
Here is a plot of k^-2 (in red), five A exp(-Bk) functions (in blue), and their sum (in green), up to k=100 on a logarithmic y axis:
Same method again, using this function
$$ A \sum_{k=1}^{\lceil n/2 \rceil} (-1)^{k-1} e^{-B(2k-1)} \sin((2k-1)\theta) $$ $$ A \cdot \text{Imag} \left\{ \sum_{k=1}^{\lceil n/2 \rceil} (-1)^{k-1} e^{-B(2k-1)} e^{(2k-1)i\theta} \right\} $$ $$ A \cdot \text{Imag} \left\{ \sum_{k=1}^{\lceil n/2 \rceil} (-1)^{k-1} e^{(-B+i\theta)^{(2k-1)}} \right\} $$ $$ \text{Let} \quad m=\{n \quad \text{if n is odd or} \quad n-1 \quad \text{if n is even}\} $$ $$ \text{Let} \quad t = (-1)^{(m-1)/2} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - t \cdot e^{(m+2)(-B+i\theta)}} {1 + e^{2(-B+i\theta)}} \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ ( e^{-B+i\theta} - t \cdot e^{(m+2)(-B+i\theta)}) \cdot (1+e^{2(-B-i\theta)}) } {(1 + e^{2(-B+i\theta)})(1+e^{2(-B-i\theta)})} \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B+i\theta} - t \cdot e^{(m+2)(-B+i\theta)} + e^{2(-B-i\theta)} e^{(-B+i\theta)} - t \cdot e^{2(-B-i\theta)} e^{(m+2)(-B+i\theta)} } { 1 + e^{-2B} (e^{2i\theta} + e^{-2i\theta}) + e^{-4B} } \right \} $$ $$ A \cdot \text{Imag} \left \{ \frac{ e^{-B} e^{i\theta} - t \cdot e^{-(m+2)B} e^{(m+2)i\theta} + e^{-3B} e^{-i\theta} - t \cdot e^{-(m+4)B} e^{(m)i\theta} } { 1 + e^{-2B} (e^{2i\theta} + e^{-2i\theta}) + e^{-4B} } \right \} $$ $$ A \cdot \frac{ e^{-B} \sin(\theta) - t \cdot e^{-(m+2)B} \sin((m+2)\theta) - e^{-3B} \sin(\theta) - t \cdot e^{-(m+4)B} \sin(m\theta) } { 1 + 2e^{-2B} cos(2\theta) + e^{-4B} } $$ $$ A \cdot \frac{ (e^{-B}- e^{-3B}) \sin(\theta) - (-1)^{(m-1)/2} \cdot (e^{-(m+2)B} \sin((m+2)\theta) + e^{-(m+4)B} \sin(m\theta)) } { 1 + 2e^{-2B} cos(2\theta) + e^{-4B} } $$Putting that back into our square approximation, we get (where m is odd)
$$ \text{triangle}_m(\theta) = $$ $$ \frac{2}{\pi} \sum_{k=1}^{5} A_k \cdot \frac{ (e^{-B_k}- e^{-3B_k}) \sin(\theta) - (-1)^{(m-1)/2} \cdot (e^{-(m+2)B_k} \sin((m+2)\theta) + e^{-(m+4)B_k} \sin(m\theta)) } { 1 + 2e^{-2B_k} cos(2\theta) + e^{-4B_k} } $$Here is a plot of that approximation at n=15
For comparison, here is the exact triangle wave at n=15