Experimenting with Chebyshev polynomials to do additive synthesis based on a single sine and cosine value pair.
Rough estimate 5 floating point arithmetic operations for each additional next harmonic, of which 3 are for the polynomial recurrence and 2 are for gain and accumulation. Compared to 3 flops and 1 `sin()` call for the obvious way. The `sin()` call way can do sparse harmonics much more easily though.
Oh! Cosines diverge as sum(1/k) diverges. Alternating signs just pushes the problem to the opposite phase (+/-PI).
Implemented bandlimited PWM via difference of two bandlimited saws (constructed by additive synthesis):
```
vec2 saw(float t)
{
float s = sin(t);
float c = cos(t);
float su = 1.0;
float sv = 2.0 * c;
float cu = 1.0;
float cv = c;
float sgn = 1.0;
int k = 1;
float sum = 0.0;
float dsum = 0.0;
while (k <= 16)
{
float term = sgn * su / float(k++);
float dterm = sgn * cu;
sum += term;
dsum += dterm;
sgn = -sgn;
float sw = 2.0 * c * sv - su; su = sv; sv = sw;
float cw = 2.0 * c * cv - cu; cu = cv; cv = cw;
}
return vec2(s * sum, dsum);
}
```
(This computes the derivative too for waveform plotting purposes.)
If you need a lot of harmonics, it can be done in parallel (SIMD? FPGA?) by using the relations:
$$
T_{2n} = T_n^2 - 1
T_{2n+1} = T_{n+1} T_n - x
$$
For the second kind, I think you can use these ones:
$$
U_{2n-1} = 2 T_n U_{n-1}
U_{2n} = T_{2n} + x U_{2n-1}
$$
(hope I didn't make any mistakes working that out)
Cosine harmonics of the form 2k+1 (eg for square wave) can be generated with the same efficiency as cosine harmonics of the form k (eg for sawtooth wave):
$$ T_{k+2} = 2 T_1 T_{k+1} - T_{k} $$
vs
$$ T_{2(k+2)+1} = 2 T_2 T_{2(k+1)+1} - T_{2k+1} $$
(the waveforms are defined with sines, but maybe cosines wouldn't sound too different)