## The Digital Waveguide Oscillator

In this section, adapted from [460], a digital
sinusoidal oscillator derived from digital waveguide theory is
described which has good properties for VLSI implementation. Its main
features are *no wavetable* and a computational complexity of
only *one multiply* per sample when amplitude and frequency are
constant. Three additions are required per sample. A piecewise
exponential amplitude envelope is available for the cost of a second
multiplication per sample, which need not be as expensive as the
tuning multiply. In the presence of frequency modulation (FM), the
amplitude coefficient can be varied to exactly cancel amplitude
modulation (AM) caused by changing the frequency of oscillation.

### Additive Synthesis

One of the very first computer music techniques introduced was *additive
synthesis* [379]. It is based on Fourier's theorem which
states that any sound can be constructed from elementary sinusoids, such as
are approximately produced by carefully struck tuning forks. Additive
synthesis attempts to apply this theorem to the synthesis of sound by
employing large banks of sinusoidal oscillators, each having independent
amplitude and frequency controls. Many analysis methods, *e.g.*, the phase
vocoder, have been developed to support additive synthesis. A summary is
given in [424].

While additive synthesis is very powerful and general, it has been held back from widespread usage due to its computational expense. For example, on a single DSP56001 digital signal-processing chip, clocked at 33 MHz, only about sinusoidal partials can be synthesized in real time using non-interpolated, table-lookup oscillators. Interpolated table-lookup oscillators are much more expensive, and when all the bells and whistles are added, and system overhead is accounted for, only around fully general, high-quality partials are sustainable at KHz on a DSP56001 (based on analysis of implementations provided by the NeXT Music Kit).

At CD-quality sampling rates, the note A1 on the piano requires sinusoidal partials, and at least the low-frequency partials should use interpolated lookups. Assuming a worst-case average of partials per voice, providing 32-voice polyphony requires partials, or around DSP chips, assuming we can pack an average of partials into each DSP. A more reasonable complement of DSP chips would provide only -voice polyphony which is simply not enough for a piano synthesis. However, since DSP chips are getting faster and cheaper, DSP-based additive synthesis looks viable in the future.

The cost of additive synthesis can be greatly reduced by making special
purpose VLSI optimized for sinusoidal synthesis. In a VLSI environment,
major bottlenecks are *wavetables* and *multiplications*. Even
if a single sinusoidal wavetable is shared, it must be accessed
sequentially, inhibiting parallelism. The wavetable can be eliminated
entirely if *recursive algorithms* are used to synthesize sinusoids
directly.

### Digital Sinusoid Generators

In [168], three techniques were examined for
generating sinusoids digitally by means of recursive
algorithms.^{C.12} The recursions can be
interpreted as implementations of second-order digital resonators in
which the damping is set to zero. The three methods considered were
(1) the
*2D rotation* (2DR), or
*complex multiply* (also called the ``coupled form''), (2) the
*modified coupled form* (MCF), or ``magic circle'' algorithm,^{C.13}which is similar to (1) but with better numerical behavior, and (3)
the direct-form, second-order, digital resonator (DFR) with its poles
set to the unit circle.

These three recursions may be defined as follows:

The digital waveguide oscillator appears to have the best overall properties yet seen for VLSI implementation. This structure, introduced in [460], may be derived from the theory of digital waveguides (see Appendix C, particularly §C.9, and [433,464]). Any second-order digital filter structure can be used as a starting point for developing a corresponding sinusoidal signal generator, so in this case we begin with the second-order waveguide filter.

### The Second-Order Waveguide Filter

The first step is to make a second-order digital filter with zero damping by abutting two unit-sample sections of waveguide medium, and terminating on the left and right with perfect reflections, as shown in Fig.C.38. The wave impedance in section is given by , where is air density, is the cross-sectional area of tube section , and is sound speed. The reflection coefficient is determined by the impedance discontinuity via . It turns out that to obtain sinusoidal oscillation, one of the terminations must provide an inverting reflection while the other is non-inverting.

At the junction between sections and , the signal is partially
transmitted and partially reflected such that energy is conserved, *i.e.*, we
have *lossless scattering*. The formula for the reflection
coefficient can be derived from the physical constraints that (1)
pressure is continuous across the junction, and (2) there is no net flow
into or out of the junction. For traveling pressure waves and
volume-velocity waves , we have
and
. The physical pressure and volume velocity are obtained by
summing the traveling-wave components.

The discrete-time simulation for the physical system of Fig.C.38 is shown
in Fig.C.39. The propagation time from the junction to a reflecting
termination and back is one sample period. The half sample delay from the
junction to the reflecting termination has been *commuted* with the
termination and combined with the half sample delay *to* the
termination. This is a special case of a ``half-rate'' waveguide filter
[433].

Since only two samples of delay are present, the digital system is at most second order, and since the coefficients are real, at most one frequency of oscillation is possible in .

The scattering junction shown in the figure is called the Kelly-Lochbaum junction in the literature on lattice and ladder digital filters [173]. While it is the most natural from a physical point of view, it requires four multiplies and two additions for its implementation.

It is well known that lossless scattering junctions can be implemented in a
variety of equivalent forms, such as the two-multiply and even one-multiply
junctions. However, most have the disadvantage of not being *normalized* in the sense that changing the reflection coefficient
changes the amplitude of oscillation. This can be understood physically by
noting that a change in implies a change in . Since the
signal power contained in a waveguide variable, say , is
, we find that modulating the reflection coefficient
corresponds to modulating the signal energy represented by the signal
sample in at least one of the two delay elements. Since energy is
proportional to amplitude squared, energy modulation implies amplitude
modulation.

The well-known normalization procedure is to replace the traveling
pressure waves by ``root-power'' pressure waves
so that signal power is just the square of a signal
sample
. When this is done, the scattering junction
transforms from the Kelly-Lochbaum or one-multiply form into the
*normalized ladder* junction in which the reflection coefficients
are again , but the forward and reverse transmission
coefficients become
. Defining
, the
transmission coefficients can be seen as
, and we arrive
essentially at the *coupled form*, or two-dimensional vector
rotation considered in [168].

An alternative normalization technique is based on the digital waveguide
*transformer* (§C.16). The purpose of a ``transformer'' is to
``step'' the force variable (pressure in our example) by some factor
without scattering and without affecting signal energy. Since traveling
signal power is proportional to pressure times velocity , it
follows that velocity must be stepped by the inverse factor to keep
power constant. This is the familiar behavior of transformers for analog
electrical circuits: voltage is stepped up by the ``turns ratio'' and
current is stepped down by the reciprocal factor. Now, since
, traveling signal power is equal to
. Therefore,
stepping up pressure through a transformer by the factor corresponds to
stepping up the wave impedance by the factor . In other words,
the transformer raises pressure and decreases volume velocity by raising
the wave impedance (narrowing the acoustic tube) like a converging cone.

If a transformer is inserted in a waveguide immediately to the left, say, of a scattering junction, it can be used to modulate the wave impedance ``seen'' to the left by the junction without having to use root-power waves in the simulation. As a result, the one-multiply junction can be used for the scattering junction, since the junction itself is not normalized. Since the transformer requires two multiplies, a total of three multiplies can effectively implement a normalized junction, where four were needed before. Finally, in just this special case, one of the transformer coefficients can be commuted with the delay element on the left and combined with the other transformer coefficient. For convenience, the coefficient on the left is commuted into the junction so it merely toggles the signs of inputs to existing summers. These transformations lead to the final form shown in Fig.C.40.

The ``tuning coefficient'' is given by
, where
is the desired oscillation frequency in Hz at sample (in the
undamped case), and is the sampling period in seconds. The
``amplitude coefficient'' is
, where growth or
decay factor per sample (
for constant
amplitude),^{C.14} and is the normalizing transformer ``turns
ratio'' given by

When amplitude and frequency are constant, there is no gradual exponential
growth or decay due to round-off error. This happens because the only
rounding is at the output of the tuning multiply, and all other
computations are exact. Therefore, quantization in the tuning coefficient
can only cause quantization in the frequency of oscillation. Note that any
one-multiply digital oscillator should have this property. In contrast,
the only other known normalized oscillator, the coupled form, *does*
exhibit exponential amplitude drift because it has *two* coefficients
and
which, after quantization, no longer
obey for most tunings.

### Application to FM Synthesis

The properties of the new oscillator appear well suited for FM applications in VLSI because of the minimized computational expense. However, in this application there are issues to be resolved regarding conversion from modulator output to carrier coefficients. Preliminary experiments indicate that FM indices less than are well behaved when the output of a modulating oscillator simply adds to the coefficient of the carrier oscillator (bypassing the exact FM formulas). Approximate amplitude normalizing coefficients have also been derived which provide a first-order approximation to the exact AM compensation at low cost. For music synthesis applications, we believe a distortion in the details of the FM instantaneous frequency trajectory and a moderate amount of incidental AM can be tolerated since they produce only second-order timbral effects in many situations.

### Digital Waveguide Resonator

Converting a second-order oscillator into a second-order filter requires merely introducing damping and defining the input and output signals. In Fig.C.40, damping is provided by the coefficient , which we will take to be a constant

*state*of the resonator. Let us denote by the output of the delay element on the left in Fig.C.40 and let be the delay-element output on the right. In general, an output signal may be formed as any linear combination of the state variables:

The foregoing modifications to the digital waveguide oscillator result
in the so-called *digital waveguide resonator* (DWR)
[304]:

where, as derived in the next section, the coefficients are given by

where denotes one desired pole (the other being at ). Note that when (undamped case). The DWR requires only two multiplies per sample. As seen earlier, when the decay time is set to (), one of the multiplies disappears, leaving only

*one*multiply per sample for sinusoidal oscillation.

Figure C.41 shows an overlay of initial impulse responses for
the three resonators discussed above. The decay factor was set to
, and the output of each multiplication was quantized to 16
bits, as were all coefficients. The three waveforms sound and look
identical. (There *are* small differences, however, which can be
seen by plotting the differences of pairs of waveforms.)

Figure C.42 shows the same impulse-response overlay but with and only 4 significant bits in the coefficients and signals. The complex multiply oscillator can be seen to decay toward zero due to coefficient quantization ( ). The MCF and DWR remain steady at their initial amplitude. All three suffer some amount of tuning perturbation.

### State-Space Analysis

We will now use state-space analysis^{C.15}[449] to determine
Equations (C.133-C.136).

or, in vector notation,

(C.137) | |||

(C.138) |

where we have introduced an input signal , which sums into the state vector via the (or ) vector . The output signal is defined as the vector times the state vector . Multiple outputs may be defined by choosing to be an matrix.

A basic fact from linear algebra is that the *determinant* of a
matrix is equal to the *product of its eigenvalues*. As a quick
check, we find that the determinant of is

When the eigenvalues of (system poles) are complex, then they must form a complex conjugate pair (since is real), and we have . Therefore, the system is

*stable*if and only if . When making a digital sinusoidal oscillator from the system impulse response, we have , and the system can be said to be ``marginally stable''. Since an undriven sinusoidal oscillator must not lose energy, and since every lossless state-space system has unit-modulus eigenvalues (consider the modal representation), we expect , which occurs for .

Note that . If we diagonalize this system to obtain , where diag, and is the matrix of eigenvectors of , then we have

If this system is to generate a real sampled sinusoid at radian frequency , the eigenvalues and must be of the form

(in either order) where is real, and denotes the sampling interval in seconds.

Thus, we can determine the frequency of oscillation (and verify that the system actually oscillates) by determining the eigenvalues of . Note that, as a prerequisite, it will also be necessary to find two linearly independent eigenvectors of (columns of ).

###

Eigenstructure

Starting with the defining equation for an eigenvector and its corresponding eigenvalue ,

We normalized the first element of to 1 since is an eigenvector whenever is. (If there is a missing solution because its first element happens to be zero, we can repeat the analysis normalizing the second element to 1 instead.)

Equation (C.141) gives us two equations in two unknowns:

Substituting the first into the second to eliminate , we get

As approaches (no damping), we obtain

They are linearly independent provided . In the undamped case (), this holds whenever . The eigenvectors are finite when . Thus, the nominal range for is the interval .

We can now use Eq.(C.142) to find the eigenvalues:

#### Damping and Tuning Parameters

The tuning and damping of the resonator impulse response are governed by the relation

To obtain a specific decay time-constant , we must have

Therefore, given a desired decay time-constant (and the sampling interval ), we may compute the damping parameter for the digital waveguide resonator as

To obtain a desired frequency of oscillation, we must solve

for , which yields

#### Eigenvalues in the Undamped Case

When , the eigenvalues reduce to

where denotes the angular advance per sample of the oscillator. Since corresponds to the range , we see that in this range can produce oscillation at any digital frequency.

For , the eigenvalues are real, corresponding to exponential growth and/or decay. (The values were excluded above in deriving Eq.(C.144).)

In summary, the coefficient in the digital waveguide oscillator () and the frequency of sinusoidal oscillation is simply

### Summary

This section introduced and analyzed the digital waveguide oscillator and resonator, as well as some related algorithms. As a recursive algorithm for digital sinusoid generation, it has excellent properties for VLSI implementation. It is like the 2D rotation (complex multiply) in that it offers instantaneous amplitude from its state and constant amplitude in the presence of frequency modulation. However, its implementation requires only two multiplies per sample instead of four. When used as a constant-frequency oscillator, it requires only one multiply per sample.

### Matlab Sinusoidal Oscillator Implementations

This section provides Matlab/Octave program listings for the sinusoidal resonator/oscillator algorithms discussed above:

- Planar 2D Rotation (2DR) (``complex multiply'')
- Modified Coupled Form (MCF) (``Magic Circle'')
- Digital Waveguide Resonator (DWR)

% Filter test program in matlab %N = 300; % Number of samples to generate N = 3000; % Number of samples to generate f = 100; % Desired oscillation frequency (Hz) fs = 8192; % Audio sampling rate (Hz) %R = .99; % Decay factor (1=never, 0=instant) R = 1; % Decay factor (1=never, 0=instant) b1 = 1; % Input gain to state variable x(n) b2 = 0; % Input gain to state variable y(n) %nd = 16; % Number of significant digits to use nd = 4; % Number of significant digits to use base = 2; % Mantissa base (two significant figures) % (see 'help chop' in Matlab) u = [1,zeros(1,N-1)]; % driving input test signal theta = 2*pi*f/fs; % resonance frequency, rad/sample % ================================================ % 2D PLANAR ROTATION (COMPLEX MULTIPLY) x1 = zeros(1,N); % Predeclare saved-state arrays y1 = zeros(1,N); x1(1) = 0; % Initial condition y1(1) = 0; % Initial condition c = chop(R*cos(theta),nd,base); % coefficient 1 s = chop(R*sin(theta),nd,base); % coefficient 2 for n=1:N-1, x1(n+1) = chop( c*x1(n) - s*y1(n) + b1*u(n), nd,base); y1(n+1) = chop( s*x1(n) + c*y1(n) + b2*u(n), nd,base); end % ================================================ % MODIFIED COUPLED FORM ("MAGIC CIRCLE") % (ref: http://ccrma.stanford.edu/~jos/wgo/ ) x2 = zeros(1,N); % Predeclare saved-state arrays y2 = zeros(1,N); x2(1) = 0.0; % Initial condition y2(1) = 0.0; % Initial condition e = chop(2*sin(theta/2),nd,base); % tuning coefficient for n=1:N-1, x2(n+1) = chop(R*(x2(n)-e*y2(n))+b1*u(n),nd,base); y2(n+1) = chop(R*(e*x2(n+1)+y2(n))+b2*u(n),nd,base); end % ================================================ % DIGITAL WAVEGUIDE RESONATOR (DWR) x3 = zeros(1,N); % Predeclare saved-state arrays y3 = zeros(1,N); x3(1) = 0; % Initial condition y3(1) = 0; % Initial condition g = R*R; % decay coefficient t = tan(theta); % toward tuning coefficient cp = sqrt(g/(g + t^2*(1+g)^2/4 + (1-g)^2/4)); % exact %cp = 1 - (t^2*(1+g)^2 + (1-g)^2)/(8*g); % good approx b1n = b1*sqrt((1-cp)/(1+cp)); % input scaling % Quantize coefficients: cp = chop(cp,nd,base); g = chop(g,nd,base); b1n = chop(b1n,nd,base); for n=1:N-1, gx3 = chop(g*x3(n), nd,base); % mpy 1 y3n = y3(n); temp = chop(cp*(gx3 + y3n), nd,base); % mpy 2 x3(n+1) = temp - y3n + b1n*u(n); y3(n+1) = gx3 + temp + b2*u(n); end % ================================================ % playandplot.m figure(4); title('Impulse Response Overlay'); ylabel('Amplitude'); xlabel('Time (samples)'); alt=1;%alt = (-1).^(0:N-1); % for visibility plot([y1',y2',(alt.*y3)']); grid('on'); legend('2DR','MCF','WGR'); title('Impulse Response Overlay'); ylabel('Amplitude'); xlabel('Time (samples)'); saveplot(sprintf('eps/tosc%d.ps',nd)); if 1 playandplot(y1,u,fs,'2D rotation',1); playandplot(y2,u,fs,'Magic circle',2); playandplot(y3,u,fs,'WGR',3); end function playandplot(y,u,fs,ttl,fnum); % sound(y,fs,16); figure(fnum); unwind_protect # protect graph state subplot(211); axis("labely"); axis("autoy"); xlabel(""); title(ttl); ylabel('Amplitude'); xlabel('Time (ms)'); t = 1000*(0:length(u)-1)/fs; timeplot(t,u,'-'); legend('input'); subplot(212); ylabel('Amplitude'); xlabel('Time (ms)'); timeplot(t,y,'-'); legend('output'); unwind_protect_cleanup # restore graph state grid("off"); axis("auto","label"); oneplot(); end_unwind_protect endfunction

Note that the `chop` utility only exists in Matlab.
In Octave, the following ``compatibility stub'' can be used:

function [y] = chop(x,nd,base) y = x;

### Faust Implementation

The function `oscw` (file `osc.lib`) implements a
digital waveguide sinusoidal oscillator in the Faust language
[453,154,170]. There is also
`oscr` implementing the 2D rotation case, and `oscs`
implementing the modified coupled form (magic circle).

**Next Section:**

Non-Cylindrical Acoustic Tubes

**Previous Section:**

Waveguide Transformers and Gyrators