Hello,
this is a very efficient way to get quasi-bandlimited sawtooth and pulse waveforms. I don't think this technique has been proposed before but I could be wrong. By no means is this idea limited to Csound, in fact I hope it will be implemented in all kinds of environments and soft synths, but since my language of choice is Csound I might just as well present it here. The following explanation is quite technical so feel free to skip right on to the examples.
The idea is to convolve an impulse train with a sin^n window (=sin(x)^n from 0 to pi) and then integrate the result. That is, the impulses are replaced with sin^n windows. The integration can be done analytically so this allows us to directly calculate the sawtooth waveform. No leaky integrators or DC blocking filters needed!
The frequency response of a sin^n window is such that the main lobe goes to zero at n/2+1 after which the sidelobes zero at intervals of 1. The main lobe gets wider as n gets larger while the side lobes are increasingly attenuated.
For n=7, I calculated the attenuation of the first four side lobe maxima relative to the main lobe maximum to be about -68, -85, -98 and -108 dB. If the zero of the main lobe is set at the Nyquist frequency this is how much the aliased harmonics are attenuated and this in addition to the normal roll-off of sawtooth wave so it's pretty much bandlimited for all practical purposes.
Mathematically I think there is some advantage in choosing an odd exponent. It simplifies the result of the integration and the expression for DC correction.
The first example is a sawtooth instrument:
sr = 44100
ksmps = 1
nchnls = 1
0dbfs = 1
instr 1
icps = p4
imainlobezero = 22050
iwindowwidth = 4.5/imainlobezero
icorrect = .5*iwindowwidth*icps
aramp phasor icps, icorrect
aindex = 1/icps/iwindowwidth*aramp
aread tablei aindex, 1, 1
aout = aread-aramp+icorrect
adeclick linseg 0, 1/icps, 1, p3-2/icps, 1, 1/icps, 0
out adeclick*aout
endin
f 1 0 65537 -9 .5 -.59814453125 90 1.5 .11962890625 90 2.5 -.02392578125 90 3.5 .00244140625 90
i 1 0 1 220
i 1 1 1 440
i 1 2 1 880
i 1 3 1 1760
e
Scaling a sin(x)^7 window to 35*pi/32*sin(pi*x)^7 and integrating gives
-1225/2048*cos(pi*x)+245/2048*cos(3*pi*x)-49/2048*cos(5*pi*x)+5/2048*cos(7*pi*x)
(x goes from 0 to 1)
The scaling normalizes the integrated window so that it goes from -1/2 at x=0 to 1/2 at x=1. This simplifies things later. The f statement in the score loads this into a lookup table.
The main lobe zero is set at 22050. The 4.5 at the expression for the window width comes from 7/2+1. Note that the main lobe zero frequency shouldn't be smaller than 4.5 times the note frequency, otherwise a new window starts before the previous one has finished. With a sampling rate of 44100 this requirement gives us a maximum allowed frequency of 22050/4.5=4900 which is sufficient for musical purposes I think.
The icorrect is an expression for both a DC correction and for an optional phase correction that sets the waveform to start from zero. A phasor at the note frequency is produced. The aindex gives us the index by which the stored lookup table is read. When the aindex value exceeds 1 the read value stays at 1/2.
Now comes the simplification that the window scaling allows: we just subtract the phasor value from the table value and add the DC correction and we get a quasi-bandlimited sawtooth, neat! The subtracted ramp acts as the integrated DC correction of the windowed impulse train.
A possible disadvantage of this method is the slightly greater than ideal roll-off of the harmonics which can be countered with a higher sampling rate or by oversampling. This also allows a higher maximum frequency if desired.
Here is a pulse instrument which works by subtracting two phase shifted sawtooths (nothing new here):
sr = 44100
ksmps = 1
nchnls = 1
0dbfs = 1
instr 1
icps = p4
imainlobezero = 22050
iwindowwidth = 4.5/imainlobezero
idutycycle = .25
aramp1 phasor icps
aramp2 phasor icps, idutycycle
ifactor = 1/icps/iwindowwidth
aindex1 = ifactor*aramp1
aindex2 = ifactor*aramp2
aread1 tablei aindex1, 1, 1
aread2 tablei aindex2, 1, 1
aout = aread1-aramp1-aread2+aramp2
adeclick linseg 0, 1/icps, 1, p3-2/icps, 1, 1/icps, 0
out adeclick*aout
endin
f 1 0 65537 -9 .5 -.59814453125 90 1.5 .11962890625 90 2.5 -.02392578125 90 3.5 .00244140625 90
i 1 0 1 220
i 1 1 1 440
i 1 2 1 880
i 1 3 1 1760
e
Note that changing the main lobe zero frequency gives us low pass filtering for free! The last example is a simulation of a resonant filter sweep:
sr = 44100
ksmps = 1
nchnls = 1
0dbfs = 1
instr 1
afeedback init 0
icps = p4
asweep expsega 1, 2, .5
amainlobezero = 4.5*icps+asweep*(22050-4.5*icps)
awindowwidth = 4.5/amainlobezero
icorrect = .5*4.5/22050
acorrect = .5*awindowwidth*icps
aramp phasor icps, icorrect
aindex = 1/icps/awindowwidth*aramp
aread tablei aindex, 1, 1
adeclick linseg 0, 1/icps, 1, p3-2/icps, 1, 1/icps, 0
aout = adeclick*(aread-aramp+acorrect-.875*afeedback)
afeedback vdelay aout, 1000/amainlobezero, 1000/icps
out adeclick*aout
endin
f 1 0 65537 -9 .5 -.59814453125 90 1.5 .11962890625 90 2.5 -.02392578125 90 3.5 .00244140625 90
i 1 0 16 55
e
Now, I doubt that these techniques can be extended quite as simply and elegantly to parabolic and triangle waves but I hope I'm wrong.
Kalle Aho
| Attachment | Size |
|---|---|
| saw.csd | 703 bytes |
| pulse.csd | 804 bytes |
| rez.csd | 851 bytes |
I realized that the program
I realized that the programs can be simplified further, we can get rid of the DC correction. Note that the lookup table now goes from -1 to 1. Here is the program for sawtooth:
sr = 44100
ksmps = 1
nchnls = 1
0dbfs = 1
instr 1
icps = p4
imainlobezero = 22050
aramp phasor icps
ascaledramp = 2*aramp-1
aindex = .5*imainlobezero/4.5/icps*ascaledramp
aread tablei aindex, 1, 1, .5
aout = aread-ascaledramp
adeclick linseg 0, 1/icps, 1, p3-2/icps, 1, 1/icps, 0
out adeclick*aout
endin
f 1 0 65537 -9 .5 -1.1962890625 90 1.5 .2392578125 90 2.5 -.0478515625 90 3.5 .0048828125 90
i 1 0 1 220
i 1 1 1 440
i 1 2 1 880
i 1 3 1 1760
e
Since I've received some
Since I've received some questions about the technique I'd like to present the idea more clearly now. There
actually are no sine windows, impulse trains or integrations in the program itself, they are just the
mathematical scaffolding behind the technique. You can ignore those and just think of the sawtooth generation
like this:
the lookup table contains a sigmoid that goes from -1 at -1 to 1 at 1.
http://en.wikipedia.org/wiki/Sigmoid_function
The particular sigmoid I use is chosen because it has predictable spectral properties, the frequency where
the magnitude spectrum goes to zero is easy to calculate. You could use something like truncated tanh() as
Victor Lazzarini pointed out to me in the Csound list.
Now think of the phasor at the note frequency as going from -1 to 1. You get the index that reads the lookup
table by multiplying the phasor with a>=1. Thus, the table reading index goes from -a to a. When the index is
between -1 and 1 the lookup table is read normally but when the index is <-1 the value stays at -1 and when
it is >1 the value stays at 1.
So when a>=1 you always get a sigmoid going from -1 to 1 repeatedly at the note frequency. Now the program
subtracts the phasor from that. This produces the sawtooth. Varying the variable a gives you a different
number of harmonics. It's quite simple really.
* * *
But now for something remarkable!
We can think of this technique more generally. If the sigmoid is an odd function centered at zero (like tanh
for example) that is defined from -inf to inf and increases monotonically, when x goes from -1 to 1 we get an
approximate sawtooth with
f(x)=x-sigmoid(a*x)/sigmoid(a)
now let's try it with tanh:
f(x)=x-tanh(a*x)/tanh(a)
when x goes from -1 to 1, we get an approximate sawtooth that has roughly 2*a harmonics at least when a is
about >=3.5
We don't have to use a lookup table, this sawtooth can be computed directly! I'm not sure how fast it is
though.
Now this can be integrated to get
f(x)=1/2*x^2-ln(cosh(a*x))/(a*tanh(a))
which gives a parabolic wave when x goes from -1 to 1. Two phase shifted parabolic waves can be subtracted to
get variable duty cycle triangle waves.
Kalle Aho