Virtual Analog Models of the Lockhart and Serge Wavefolders †

: Wavefolders are a particular class of nonlinear waveshaping circuits, and a staple of the “West Coast” tradition of analog sound synthesis. In this paper, we present analyses of two popular wavefolding circuits—the Lockhart and Serge wavefolders—and show that they achieve a very similar audio effect. We digitally model the input–output relationship of both circuits using the Lambert-W function, and examine their time-and frequency-domain behavior. To ameliorate the issue of aliasing distortion introduced by the nonlinear nature of wavefolding, we propose the use of the ﬁrst-order antiderivative method. This method allows us to implement the proposed digital models in real-time without having to resort to high oversampling factors. The practical synthesis usage of both circuits is discussed by considering the case of multiple wavefolder stages arranged in series.


Introduction
Nonlinear waveshaping is a technique used in sound synthesis to generate complex harmonic spectra.It consists of processing a signal with low harmonic content (typically a sinusoid) using a nonlinear mapping function designed to introduce harmonic overtones to the output signal [1].The first documented use of waveshaping in the digital domain can be traced back to 1969, when Jean-Claude Risset emulated the sound of a clarinet by distorting a sinusoid with a clipping function [2].Waveshaping techniques were extensively researched within the context of computer music in the 1970s, with several authors exploring the use of Chebyshev polynomials in particular, as an accurate and computationally cheap alternative to additive synthesis [1,[3][4][5].The underlying principles behind waveshaping synthesis are closely related to other well-known synthesis techniques, such as frequency modulation (FM) and phase distortion (PD) synthesis [6,7].These two techniques rely on distorting the frequency and phase, respectively, of sinusoidal oscillators.Recent research on the topic of distortion-based synthesis has explored the use of logic operators in lieu of traditional polynomial waveshapers [8], and proposed extensions to both FM and PD synthesis [9,10].
The use of waveshaping in the analog domain began in the 1950s, when guitar players started deliberately overdriving their tube amplifiers to alter the timbre of their instrument [11].In 1961, Gibson released the "Maestro FZ-1 Fuzz Tone", the first commercially available fuzz distortion pedal, which exploited the saturating behavior of transistors to introduce harmonic distortion [12].Most guitar distortion pedals, including popular designs such as the Ibanez Tube Screamer and Electro-Harmonix Big Muff Pi, operate under this same basic principle [13,14].
In analog synthesizers, the use of distortion-based methods is one of the cornerstones of "West Coast" synthesis, a paradigm pioneered by California-native Don Buchla during the 1960s.Buchla's instruments focused on timbre manipulation at oscillator level by employing a variety of techniques such as nonlinear waveshaping, oscillator synchronization and pitch modulation [15][16][17].This approach to sound synthesis contrasts that of traditional subtractive synthesis, where timbre is typically controlled by filtering harmonically-rich oscillator waveforms, like sawtooth and square waves, using resonant filters [18].In recent years, West Coast synthesis has become increasingly popular, with contemporary manufacturers such as Make Noise and Verbos Electronics releasing their own takes on classic Buchla circuits.
This study presents virtual analog (VA) models for two analog synthesizer circuits: the Lockhart wavefolder and the wavefolder used in the middle section of the Serge Wave Multipliers.Wavefolding is a type of nonlinear waveshaping common in West Coast synthesis where portions of the input signal that exceed certain threshold are inverted or "folded back", hence the name of the effect.The two circuits considered in this study were chosen because of the strong similarities between their general behavior.In a similar way to guitar distortion pedals, both wavefolders exploit the saturating behavior of semiconductor p-n junctions (i.e., transistors/diodes) to implement a folding function.
Wavefolders are amongst the most emblematic building blocks of West Coast synthesis.In spite of that, they have been mostly overlooked by both VA and digital waveshaping research.We have recently begun to fill this research gap in [17], which presents a VA model of the wavefolder circuit in the seminal Buchla 259 module.Previous work on circuit-based VA modeling has researched the filters found in vintage synthesizers such as those produced by Moog [19][20][21][22], Electronic Music Studios (EMS) [23,24], Korg [25,26] and Buchla [16].Extensive work has also been done on modeling guitar distortion pedals [13,27], tube amplifiers [11,28,29], modulation effects [30][31][32][33] and the Roland TR-808 drum machine [34,35].Measurement-based VA modeling, commonly known as "black-box modeling", has also been thoroughly studied within the context of guitar amplifiers and pedals [36][37][38].This approach is particularly useful when the original circuit schematics are not available.
A major challenge in VA modeling of nonlinear circuits, and digital waveshaping in general, is aliasing suppression.Early research on waveshaping synthesis addressed this issue by using low-order polynomial transfer functions, which not only allowed full parametric control of the produced spectrum but also ensured that the output waveform was bandlimited [4].In VA modeling, high oversampling factors are usually necessary to prevent harmonics introduced by nonlinearities from reflecting into the baseband as aliases [13].Oversampling increases the computational requirements of the model, by introducing additional filtering stages and scaling the number operations required to compute each output sample.For VA models that require evaluating transcendental functions, as is the case with the proposed Lockhart and Serge models, these added costs could compromise the integration of the system within a larger, real-time computer music system.
A sizable portion of VA research has concentrated on designing efficient algorithms to generate alias-free geometric waveforms like those used in analog subtractive synthesizers, the so-called classic analog waveforms.Well-known techniques include the bandlimited impulse train (BLIT) family of methods, which involves the use of bandlimited basis functions and their integrated forms [39][40][41], and the use of differentiated polynomial waveforms (DPW) [42][43][44].Moreover, Välimäki and Franck have applied the antialiasing principle behind the DPW algorithm to tackle aliasing in wavetable oscillators [45].Recent work on antialiasing techniques has extended the use of the bandlimited ramp (BLAMP) method, originally proposed to antialias triangular oscillators in [25,41], to special cases of linear piecewise nonlinearities such as signal rectification, and inverse/hard clipping [17,46,47].
In this work, we propose the use of the antiderivative antialiasing method introduced in [48,49].This approach can be used to reduce the aliasing caused by arbitrary nonlinear waveshaping functions and is applicable to the proposed wavefolder models.In its first-order form, the method can be derived by analytically convolving a linear continuous-time representation of the input signal with a rectangular lowpass kernel [48].As shown in this work, the use of the antiderivative method reduces the oversampling requirements of the proposed wavefolder models.
A VA model of the Lockhart wavefolder was originally presented in [50].This paper extends that work by introducing a second wavefolding circuit and studying the similarities between both systems.Additionally, we present a different treatment of the required Lambert-W function and an extended evaluation of the proposed antialiasing method in terms of computational costs.
This paper is organized as follows.Sections 2 and 3 describe the model derivation of the Lockhart and Serge wavefolders, respectively.Time-domain simulations of the circuits are also presented in these two sections.Section 4 deals with two implications of VA wavefolding in the digital domain, namely aliasing suppression and evaluation of the Lambert-W function.Section 5 presents frequency-domain results of the Lockhart and Serge wavefolders, as well as an evaluation of the proposed antialiasing method in terms of perceived sound quality and computational costs.Section 6 discusses the practical synthesis usage of both circuits and compares the behavior of the middle Serge Wave Multiplier with a recommended four-stage topology built around the Lockhart wavefolder.Concluding remarks and perspectives appear in Section 7.

The Lockhart Wavefolder
Figure 1a shows a simplified circuit diagram of the Lockhart wavefolder.This circuit was designed in 1973 by R. Lockhart Jr., who intended it to be used as a general-purpose frequency tripler [51].Following its publication in Bernie Hutchin's Electronotes [52], Lockhart's design was repurposed as a wavefolder by Ken Stone, who realized its potential as a simple yet interesting waveshaper [53,54].The Lockhart wavefolder has become ubiquitous in the music synthesizer do-it-yourself (DIY) community.For example, it is the core processor in Yves Usson's "Metalizer" module [55].
(a) (a) The main modifications made by Stone to Lockhart's original design were the addition of an input potentiometer to attenuate the amplitude of the input waveform, and an inverting amplifier at the output of the circuit [54].For the sake of simplifying the analysis, these are not shown in Figure 1a.The inverting stage at the output is reintroduced in Section 2.3.In our treatment of the circuit, we introduce the load resistance R L as an additional parameter which can be used to further control the timbre of the folded waveform.

Circuit Analysis
The Lockhart wavefolder consists of an NPN and a PNP bipolar junction transistors connected at their base and collector terminals.In order to model the large-signal behavior of the circuit, we replace transistors Q 1 and Q 2 with their corresponding Ebers-Moll large-signal models [23,24].Figure 1b shows the large-signal equivalent circuit of the Lockhart wavefolder.We use a double subscript notation to distinguish the voltages and currents in transistor Q 1 from those in Q 2 .For example, I ED,2 denotes the current through the base-emitter diode in Q 2 .Component values for the circuit have been compiled in Table 1.
We begin the analysis of the circuit by assuming that the supply voltages are equal but opposite in sign (i.e., V CC = −V EE ), and that |V in | V CC .This assumption means that the base-emitter junctions of both Q 1 and Q 2 will always be forward-biased and their voltage drops will remain approximately constant for all expected values of V in .In Ken Stone's version of the circuit, V in is assumed to be bounded between approximately ±1.2 V [53].Applying KVL around both input-emitter loops gives us where I E,1 and I E,2 are the emitter currents, and V BE,1 and V BE,2 are the voltages across the base-emitter junctions of Q 1 and Q 2 , respectively.Solving Equations ( 1) and (2) for I E,1 and I E,2 gives us: Next, we apply KCL at the collector node, which gives us where If we then assume that the contributions of the reverse currents α R I CD 1 and α R I CD 2 to the total currents associated with the emitter nodes are negligible (i.e., α R ≈ 0), we can establish that Assuming α F = 1, as suggested in [23], and inserting Equations ( 8) and ( 9) into Equations ( 6) and (7), respectively, yields a new expression for the total output current of the circuit: We then combine Equations ( 3) and (4) to derive an expression for the difference between emitter currents: Since V CC + V EE = 0, and voltage drops V BE,1 and V BE,2 are assumed to be constant and equal, their contribution to this expression disappears.Therefore, we can simplify this result as: Substituting Equation (12) into Equation ( 10) produces an expression for the total output current I out in terms of the input voltage and the currents through the collector diodes: The current-voltage (I-V) relationship of diodes can be modeled using Shockley's ideal diode equation, defined as where I D is the current through the diode, I s is the reverse bias saturation current, V D is the voltage across the diode, V T is thermal voltage and η is the ideality factor of the diode [56].For the p-n junctions inside transistors we can assume a reverse saturation current value I s = 10 −17 A and an ideality factor η = 1.A thermal voltage value V T = 25.864mV is used throughout this study.Applying Shockley's diode equation to the collector diodes and substituting into Equation (13) gives us: Next, we use KVL to derive expressions for V CD,1 and V CD,2 in terms of V in and V out : Now, the collector diodes in the large-signal model are antiparallel.Therefore, we can make the further assumption that only one of them will conduct at a time depending on the polarity of V in .A similar treatment is presented in [57] for the case of diode pairs in guitar distortion circuits.This means that By combining these new assumptions with Equations ( 15)-( 17) we arrive at the piecewise expression where λ = sgn (V in ) and sgn () Equation ( 18) can be further simplified if we consider that the independent constant factor λI s that results from its expansion will be very small (±10 −17 A) and can therefore be neglected.This gives us: Finally, we multiply both sides of this expression by R L to derive an input-output voltage relationship for the Lockhart wavefolder:

Explicit Formulation
Equation ( 21) describes an implicit relationship between the input and output voltages of the circuit; it cannot be solved algebraically.Instead, a closed-form solution for V out can be derived with the help of the Lambert-W function.The use of the Lambert-W function W() has been previously researched within the context of VA modeling.Several authors have used it to solve the implicit I-V relationship of diodes [56,58,59].Parker and D'Angelo used W() to model the Buchla Lowpass-Gate, a synthesizer circuit that employs a resistive opto-isolator (also known as a vactrol) in its control path [16].Strictly speaking, W() is multivalued; however, in this work, we only utilize the upper branch of the function.This branch is known as W 0 () in the literature [56,60].
The Lambert-W function is used to solve equations of the form which have the explicit solution where A, B, C and D ∈ R [58].Equation ( 21) can be arranged in the form described by Equation ( 22) by first rewriting it as and dividing both sides by exp (−λV out /ηV T ), which gives us Solving for V out as defined in Equation ( 23) yields an explicit model for the Lockhart wavefolder where Table 2 summarizes all parameter values for the proposed Lockhart wavefolder model.

Model Discretization and Evaluation
The voltages inside the Lockhart wavefolder are time-dependent.Therefore, we can describe the continuous-time model defined by Equation ( 26) as being of the form where f () is the transfer function of the system and t is time.In the synthesis literature, the term "transfer function" is commonly used to denote the waveshaping function [4].It should not be confused with the s-and z-domain transfer functions used in linear system analysis.
As previously mentioned, Ken Stone's circuit features an inverting stage before the output which can be modeled by inverting the polarity of the right-hand side of Equation ( 26): While including this step is not strictly necessary, we have chosen to do so, as it will facilitate the evaluation of the model.Now, given the form described by Equation ( 27), the Lockhart model can be discretized trivially by replacing all continuous-time signals with their discrete-time equivalents, i.e., where n is the sample index.
The time-domain behavior of the proposed circuit model was validated by comparing it against a reference simulation obtained using the SPICE (Simulation Program with Integrated Circuit Emphasis) software LTspice (Version IV, Linear Technology, Milpitas, CA, USA, 2016) [61].The results of this simulation are shown in Figure 2a for values of V in between −1.5 and 1.5 V. Figure 2b shows the transfer function of the proposed model implemented in MATLAB (Version R2017a, MathWorks, Natick, MA, USA, 2017) using Equation ( 28) and MATLAB's native "lambertw" function.In both cases, four different values of R L were simulated: 1, 5, 10 and 50 kΩ.From these figures, we can observe the general behavior of the Lockhart wavefolder.At low input values, the system behaves linearly, whereas for high input values the circuit inverts the slope of the driving signal.The transition between non-folded and folded portions of the signal is gradual, which responds to the characteristic soft saturating behavior of p-n junctions.The region where the transfer function folds the input signal is indicated with a blue arrow in Figure 2b for the case when R L = 50 kΩ.As shown in these Figures, increasing the value of R L sharpens the shape of the transfer function.
The curves shown in Figure 2a,b indicate a good match between the SPICE simulations and the proposed digital model.Figure 3 shows the absolute value (in mV) of the difference between both simulations.From this plot, we can observe that the difference between the curves is indeed very small, below 1 mV for all values of V in measured.These small differences are perceptually insignificant and can be attributed to the simplifications made during the analysis of the circuit and to the way in which SPICE computes currents flowing through semiconductor devices.For example, SPICE introduces a small fictitious conductance in parallel with each p-n junction in order to aid the convergence of its iterative solvers.Additionally, the SPICE diode model will account for the small reverse current that flows when the voltage across the diode is negative [61]. Figure 4 shows a time-domain view of the output of the proposed model when driven by a 500-Hz sinusoidal waveform with a peak amplitude of 1 V for two different load resistance values, R L = 10 and 50 kΩ.A sampling rate F s = 44.1 kHz was used to generate these figures, which are plotted against their corresponding SPICE simulations.From these results, we can once again observe the effect of wavefolding and the impact of R L on the overall shape of the output.For high values of R L the transition region between folded and non-folded values becomes very small and the resulting waveform is almost discontinuous (see Figure 4b).In the frequency domain, this will translate to higher harmonic content, similar to that of a square wave oscillator.A more detailed frequency-domain analysis of the Lockhart circuit is presented in Section 5 of this study.

The Serge Middle Wave Multiplier
The second circuit considered in this study is the middle section of the Serge Wave Multipliers (often abbreviated as the Serge VCM).The Serge VCM is a synthesizer module designed in 1977 by West Coast designer Serge Tcherepnin, founder of Serge Modular Music Systems.It offered three separate and independent analog sound processors, namely the "top", "middle" and "bottom" sections.As described in an original 1980 Serge product catalog, "The middle section generates a sweep of the odd harmonics (1, 3, 5, 7, 9, 11 and 13th) when a triangle wave is applied to its input...This module can be used to explore timbral areas beyond the range of ring modulation because there are more varied harmonics than the sum and difference tones" [62].
The middle Serge VCM is essentially a waveshaping circuit consisting of six identical wavefolding stages arranged in series.An amplifier at the input of the circuit is used to modulate the gain of the input waveform and control the amount of folds introduced [63].In this section we focus on the analysis of a single folding stage.The transfer function and frequency-domain behavior of the complete system are presented in Section 6. Figure 5 shows the schematic of a single wavefolding stage in the circuit.Component information is given in Table 3.To derive the transfer function for the Serge wavefolding circuit, we first assume ideal op-amp behavior and derive an expression for V out in terms of V in and V x , the voltage at the non-inverting input of the amplifier.This gives us: Since in this case R 3 = R 2 , we can further simplify this result as: Next, we derive an expression for V x by considering the subcircuit shown in Figure 6, which is essentially a diode pair similar to those found in guitar distortion circuits [13,57,58].
Equivalent view of the diode saturator at the non-inverting input of the op-amp in Figure 5.
Applying KVL around the outer loop of the circuit yields the relation where I is the current through resistor R 1 .Then, we apply KCL at the output node of the circuit, which gives us Combining Equation (32) with Equation ( 33) and applying Shockley's diode equation gives us As before, we assume the diodes will not conduct simultaneously and arrive at the piecewise relationship where once again λ = sgn(V in ).To further simplify this expression we neglect the constant factor λR 1 I s that results from its expansion.This gives us: Next, we rearrange this equation in the Lambert-W form described by Equation ( 22) by dividing both sides by exp (λV x /ηV T ).This yields which can be solved for V x as: As a final step, we insert Equation (38) into Equation (31) to derive a complete expression for the transfer function of a single wavefolding stage in the Serge middle VCM: Figure 7 shows the transfer function of the circuit, evaluated in MATLAB for values of V in between −1.5 and 1.5 V.As before, the model was discretized trivially and is presented against its corresponding SPICE simulation.Parameter values used in this simulation are compiled in Table 4.The value of parameters I s and η for the 1N4148 diode were matched to those of its corresponding SPICE model [61].Figure 7b shows the absolute difference between both simulations.These results indicate a good match between the models, as the maximum difference was once again found to be below 1 mV.Finally, Figure 8 shows the output of the Serge wavefolder for a 500-Hz sinusoidal input.As expected, the circuit behaves as a wavefolder, folding portions of the input waveform whose absolute value exceeds approximately 0.3 V.This behavior is similar to that of the Lockhart wavefolder (cf. Figure 4a).

Model Equivalence
Equation ( 39) shares a close resemblance with Equation ( 28), the proposed Lockhart wavefolder model.In fact, both expressions have the same form, which consists of the difference between a portion of the input signal and an input-dependent nonlinear element.In the case of the Lockhart wavefolder, when the R L = R/2 Equation (28) simplifies to which is remarkably close to Equation (39), with the only difference being the missing factor of two outside the Lambert-W function.This factor accounts for the difference between physical parameters I s and η in each circuit.Figure 9 shows a comparison of the transfer functions for the Lockhart (R L = 7.5 kΩ) and Serge wavefolders implemented using the parameter values in Tables 2 and 4, respectively.From this figure, we can observe that the only significant difference between both transfer functions is in their sharpness at the folding points.This means the Lockhart wavefolder will introduce sharper folds which will translate into brighter sounds at the output.From this analysis, it is clear that both circuits result in a similar audio effect, even though they are produced using different architectures.

Wavefolding in the Digital Domain
In the previous sections, the time-domain behavior of the Lockhart and Serge wavefolder models was examined via trivial discretization of their characteristic transfer functions.In this section, we move on to consider two implications of virtual analog wavefolding: evaluation of the Lambert-W function and aliasing.

Evaluating the Lambert-W Function
A particular challenge of using the Lambert-W function in VA modeling, where real-time operation is paramount, is that of computational efficiency.Optimizing the evaluation of W() is an active research topic (see, e.g., [64]).For the case of guitar distortion circuits, Paiva et al. proposed the use of a simplified iterative method which relies on a lookup table for its initial guess [57].In this work, we propose approximating the value of W() directly using Fritsch's iteration, as suggested in [60].In order to compute w m , an approximation to W(x), where x ∈ R >0 , we iterate over where and m = 0, 1, 2, . . ., M − 1.The value w 0 is an initial guess and M is the number of iterations required for ε m to approximate zero within machine-size floating point precision.The special case W(0) = 0 is defined separately.The efficiency of Fritsch's iteration will depend on the choice of initial guess.As explained in [60], an initial guess within 10 −4 of the solution will yield an approximation to W() accurate to within 10 −16 in just one iteration.Figure 10 shows the approximate times required to compute W(x) for a set of values of x between 10 −24 and 10 300 using Fritsch's iteration and the previously-proposed Halley's method [50].This range was chosen as it covers all values of interest.For instance, when V in = 5 V, the argument of W() in the Serge wavefolder model will be approximately 1.52 × 10 44 .All times were computed by averaging the result of 30 iterations implemented under identical circumstances.A piecewise approximation was used to compute the initial guess, as described in [60].From this plot, we can observe that Fritsch's iteration outperforms Halley's method by up to approximately 11 times.A MATLAB implementation of the Lambert-W function used to perform these measurements can be found in the accompanying website for this article.
Halley Fritsch Figure 10.Averaged processing times required to compute W(x) using Halley's method and Fritsch's iteration.

Aliasing Considerations
As discussed in Section 1, nonlinear waveshaping in the digital domain is susceptible to aliasing distortion due to its frequency-expanding nature.Wavefolding is no exception to this problem.As an arbitrary input waveform is folded, new harmonic overtones will be added to the frequency spectrum.Harmonics at frequencies exceeding half of the sampling rate, or the Nyquist limit, will be reflected into the audio baseband as aliases.Aliasing is known to cause unpleasant artifacts-such as beating and inharmonicity-that cannot be tolerated in a music computing scenario.Oversampling is commonly employed to mitigate this issue; however, this approach increases the computational requirements of the system by introducing additional operations.
We propose the use of the first-order antialiasing method presented in [48,49].This method is designed to reduce aliasing caused by memoryless waveshaping functions with the form described by Equation (29).The antialiased output of the waveshaping function is defined as where F() is the antiderivative of f (), the original transfer function.For the case of the Lockhart wavefolder defined by Equation ( 28), the integrated transfer function is defined as where and α, β and ∆ remain as before.This result showcases an advantageous property of the Lambert-W function W(); its antiderivative is defined in terms of W() itself.Therefore, computing F() does not pose a major increase in computational costs with respect to evaluating simply f ().For the case of the Serge wavefolder defined by Equation ( 39), the required antiderivative is given by where , Equation ( 45) can become ill-conditioned.This is avoided by defining the special case | is smaller than a predetermined threshold [48].This special case simply bypasses the antialiased form while compensating for the half-sample latency of the method.

Results
This section examines the frequency-domain behavior of the Lockhart and Serge wavefolders and their proposed antialiased forms.Next, we evaluate the computational costs of the antiderivative method with respect to oversampling for the case of the Lockhart model.

Frequency-Domain Behavior
The spectrogram in Figure 11 shows the effect of increasing the value of R L in the Lockhart wavefolder model for a constant 150-Hz sinusoidal input.As expected, the level of harmonic distortion introduced by the circuit is proportional to the value of this resistance.Therefore, this parameter can be used for additional timbral control.It should be noted that due to the antisymmetric nature of the folding function, the system introduces odd harmonics only for input signals centered around zero.Since the level of harmonics introduced by the Lockhart wavefolder depends on the choice of R L , we consider the highest recommended case R L = 50 kΩ as a worst-case scenario in terms of aliasing distortion.Figure 12 shows the spectrograms of a linear sweep from 20 Hz-5 kHz with peak amplitude of 1 V processed by the proposed Lockhart and Serge wavefolder models.This frequency range was chosen as it covers all fundamental frequencies of musical interest.A sample rate F s = 1 MHz was used to generate these figures in order to simulate an ideal alias-free continuous-time behavior.These results will be used as a reference when evaluating the performance of the proposed antialiased forms.From these spectrograms we can observe how the Lockhart wavefolder is capable of generating brighter sounds.This perceptual attribute can be varied by changing the value of R L . Figure 13a shows the result of processing the same linear sweep using the trivial (i.e., non-antialiased) Lockhart model at a standard audio rate of F s = 44.1 kHz.When compared to Figure 12a, we can clearly observe the high levels of aliasing distortion introduced by the model.This is somewhat ameliorated in Figure 13b, where the sweep has been processed at a sample rate of F s = 88.2 kHz (i.e., two-times oversampling).Figure 13c,d shows the result of processing the sweep using the proposed antialiasing method at audio rate and with two-times oversampling, respectively.As shown in these spectrograms, there is a significant reduction in aliasing, particularly below the fundamental frequency.This behavior is advantageous in music applications because at low frequencies the audibility of aliasing distortion is only limited by the hearing threshold.On the other end of the spectrum, the masking effects of harmonics will help suppress the audible effects of high-frequency aliases [65].
The spectrograms in Figure 14a,b show the outcome of processing the 1 V linear sweep with the proposed Serge wavefolder model at audio rate and with two-times oversampling, respectively.When compared with Figure 13a,b it is evident that the Serge wavefolder model is less susceptible to aliasing distortion.This can be attributed to the fact that its transfer function is not as sharp as that of the Lockhart, particularly when R L = 50 kΩ.Figure 14c,d shows the result of processing the linear sweep using the antiderivative method at audio rate and with oversampling by two.In this case, operating at audio rate yields very effective results as there are very few visible aliases left below the fundamental.When combined with oversampling by two the antiderivative method produces a nearly alias-free spectrum for the measured frequency range.The performance of the proposed antialiasing method was further evaluated by computing the A-weighted noise-to-mask ratio (ANMR) for a set of sinusoidal input signals processed by both wavefolder models.The ANMR has been previously researched as a perceptually-informed measure to evaluate the audibility of aliasing distortion [41,65].The algorithm computes the power ratio in decibels between the wanted harmonics and aliased components, but takes into account the masking effects of the former.An A-weighting filter is applied to all signals prior to the evaluation of the ANMR in order to account for the frequency-dependent sensitivity of hearing for low-level sounds.Signals with an ANMR value below −10 dB are considered to be completely free from perceivable aliasing.A detailed account of this method can be found in [65].
Figure 15a compares the measured ANMRs for a set of sinusoidal inputs with fundamental frequencies between 1 and 5 kHz processed by the Lockhart model at different sampling rates.The ideal alias-free signals required to compute these values were synthesized using Fourier analysis and additive synthesis, as suggested in [41].All signals were downsampled back to audio rate (i.e., 44.1 kHz) prior to evaluation.A dashed horizontal line has been used to indicate the −10 dB hearing threshold for aliasing distortion.In Figure 15a we can observe the significant increase in signal quality obtained by the proposed antialiasing method when applied to the Lockhart wavefolder, even when operating at audio rate.Moreover, these measurements show that the performance of the proposed method, when combined with two-times oversampling, is on par with oversampling by a factor of 8.For all fundamental frequencies below approximately 4.2 kHz, the ANMR lies below the −10 dB line.This range can be regarded as sufficient for musical applications if we consider that the highest fundamental frequency on a standard grand piano is 4186.01Hz (MIDI note C8). Figure 15b shows the measured ANMRs for the Serge wavefolder.When implemented at audio rate, the output is free from perceivable aliasing for fundamental frequencies up to approximately 2 kHz.These measurements go in accordance with the Spectrogram in Figure 14a, which shows aliasing is significantly more evident above this frequency.The use of the antiderivative method yields results comparable to those of oversampling by a factor of two, with all measured fundamental frequencies below approximately 4.6 kHz lying below the −10 dB aliasing threshold.Overall, these results indicate the proposed Serge wavefolder model can operate at audio rate with the help of the antiderivative method, therefore avoiding the need for oversampling.

Computational Costs
The computational costs of the antialiased Lockhart wavefolder model were measured by porting the algorithms into C code using the 128-bit long double data type.Table 5 shows the computation times for a 1-s 100-Hz sinusoidal input processed using the proposed model for different peak amplitude values.These results were computed by averaging the processing times of one hundred implementations.All tests were performed under identical circumstances, using a fixed resistance value of R L = 50 kΩ, the highest recommended value.From these results we can observe that the complexity of the model does not depend on the input and that the overhead of implementing the antialiasing method is minimal.When operating at audio rate, the added computation time is approximately 1 ms for a 1-second simulation.Moreover, these time measurements show that the antialiased Lockhart model, implemented at a sample rate F s = 88.2 kHz, is approximately 3.5 times faster than oversampling by a factor of 8 (i.e., F s = 352.8kHz) and nearly twice as fast as oversampling by a factor of 4. Changing the value of R L did not affect the execution times of the algorithms.For high values of R L , long double representation is necessary to account for the large values that will result at the argument of the Lambert-W function.For smaller values, 64-bit precision will be sufficient to accommodate most input levels of interest.For instance, when R L = 7.5 kΩ the signal at the input of the proposed Lockhart wavefolder model can have a peak amplitude of up to 9 V.
The measurements in Table 5 were conducted by synthesizing all input signals at the target rates.In practical implementations, oversampling will require additional pre-and post-filtering stages that will further increase the computational costs of the system.The complexity and costs of these filtering stages will be directly proportional to the required oversampling factor.This constitutes another advantage of the proposed antialiasing method.

Practical Synthesis Usage
In practical sound synthesis applications, a single folding stage is rarely used, as the timbral variety it can produce is quite limited.Most analog designs, for example the Intellijel µFold and the aforementioned Yusynth Metalizer, employ several wavefolding stages arranged in series.The number of stages varies according to the design, but typically cascades of two to six stages are used.As mentioned in Section 3, the Serge middle VCM utilizes six identical folding stages.Figure 16 shows a simplified block diagram representation of the Serge middle VCM based on the original design [63].Blocks labeled "SWF" correspond to the proposed Serge wavefolder model.An ad hoc gain factor of four, not present in the original circuit, has been added to compensate for the scaling of the signal introduced by the cascade of wavefolders.In cascaded wavefolder structures like the one shown in Figure 16, timbral control can be achieved in two manners.The first is by adjusting the gain of the input waveform (using G S in this case).This parameter controls the amount of folds introduced, allowing the overall brightness of the sound to be varied.It can be modulated in real-time to provide articulation to a sound similar to filtering in subtractive synthesis or modulation index in FM synthesis.The second way to control timbre is by adding a dc offset to the input of the wavefolder.This breaks the aforementioned symmetry of the folding function and introduces even harmonics.When modulated by using, for example, a low-frequency oscillator, this parameter provides an effect reminiscent of pulse-width modulation.Figure 17a shows the transfer function of the Serge middle VCM model for the case of zero dc offset.This plot was generated by defining V in to have a constant value of 1 V and sweeping through values of G S between −8 and 8. Figure 17b shows the output of the Serge middle VCM when driven by a 100-Hz sinusoid with G S = 6.For simplicity, in this section, we assume the range of V in to be fixed at ±1 V; therefore, all gain modulation is done using G S only.The spectrogram in Figure 18a shows the effect of increasing G S from 0 to 6 for a 150-Hz sinusoidal input.This plot effectively depicts the rich harmonic patterns introduced by the system, which are far more complex than those introduce by traditional waveshaping methods and offer a wide timbral palette for sound synthesis.The fluctuations in energy at the fundamental and first few harmonics indicate the gain values at which each new fold is introduced.Figure 17b shows the effect of introducing a dc offset at the input of the system for a constant 200-Hz sinusoidal input.This result shows how the use of a dc offset can extend the timbral possibilities of the system even further, by introducing complex patterns consisting of both even and odd harmonics.Now, although the Lockhart wavefolder was originally designed to operate as a standalone unit, it can be adapted into a series topology with relative ease.Here, we propose using the wavefolding structure shown in Figure 19 to expand the synthesis capabilities of the Lockhart wavefolder.This design, while not based on any existing circuit, is comparable to that of the Yusynth Metalizer which also utilizes four Lockhart circuits in series [55].The following paragraphs describe the sections of this proposed topology.Its frequency-domain behavior is then examined and compared with that of the Serge middle VCM.
The blocks labeled "LWF" in Figure 19 correspond to the proposed Lockhart wavefolder model.In order for this cascade of Lockhart wavefolders to behave as expected, we need to make sure that the individual folding stages satisfy two criteria.Firstly, the individual folders must provide approximately unity gain for small input values, i.e., below the folding point, and approximately negative unity gain beyond the folding point.Secondly, each stage should start folding at the same point with respect to its individual input.We can meet these criteria with the proposed Lockhart model by selecting an appropriate value for R L and adding static gain stages before and after the folding stages.These gain blocks will also help compensate for the attenuation introduced by the folding operation.First, we choose a value of L for which the Lockhart wavefolder exhibits unity gain for small input values.Having found this resistance value, the pre-and post-gain stages can be determined by measuring the value of |V out | at exactly the folding point.The pre-gain is taken to approximately this value, and the post-gain is taken to be its inverse.In Section 3.1, it was shown that for R L = 7.5 kΩ the Lockhart wavefolder exhibits approximately unity gain below the folding point.This value leads to pre-and post-gains of approximately 1/3 and 3, respectively.
Figure 20a shows the transfer function of the proposed structure measured at the output of the post-gain block.We can observe how the folds introduced by this structure are evenly distributed, unlike those in Figure 17a.As with the Serge middle VCM, timbral control is achieved by modulating the value of G L and by adding a dc offset.The static gain blocks ensure the amplitude of the folded output is bounded between approximately ±1 V for values of G L between −10 and 10 (assuming once more that V in has a peak amplitude of 1 V). Figure 20b shows the time-domain result of processing a 100-Hz sinusoidal input with the proposed structure for G L = 10 and zero dc offset.In this particular design, additional timbral control can be achieved by modulating the value of R L .Lastly, we add two optional blocks.The first is a tanh() function after the post-gain block to model the behavior of an output buffering stage and to limit the range of the output waveform.This tanh() block can also be antialiased using the antiderivative method described by Equation (45).The antiderivative of the tanh() function is given by log(cosh()) [48].The second optional block is a static one-pole lowpass filter with a cutoff at f c = 1.3 kHz whose purpose is to act as a simple tone control.A similar static filtering stage can be found at the output of the Buchla 259 wavefolder [17].The s-domain transfer function of this filtering stage is defined as where w c = 2π f c .
Finally, we examine the time-varying behavior of the proposed structure by considering the case of a 150-Hz input sinewave with variable gain G L and dc offset.Figure 21 shows the spectrogram that results from varying G L from 0 to 15.As expected, the system introduces complex harmonic patterns similar to those shown in Figure 18a.Likewise, Figure 18b demonstrates the effect of varying G L from 0 to 10 while simultaneously increasing the level of dc offset from 0 to 5 V.This response is comparable to that of the Serge middle VCM (see Figure 18b).A real-time demo of the proposed Lockhart wavefolder topology implemented using Max/MSP and Gen is available at http://research.spa.aalto.fi/publications/papers/smc17-wavefolder.

Conclusions
In this work, we have explored the behavior of two West Coast synthesizer circuits: the Lockhart and Serge wavefolders.By means of circuit analysis, we have derived closed-form expressions for the characteristic transfer functions of both systems.These transfer functions were validated against SPICE simulations implemented using LTspice.The results obtained indicate a good match between the proposed models and their corresponding SPICE simulations.In addition to this, we observed that the behavior of both circuits is very similar, despite the fact that their designs are fundamentally different.
The issue of aliasing caused by wavefolding in the digital domain was treated by incorporating the first-order antiderivative method.Within the context of the Lockhart wavefolder, it was shown that the proposed antialiased model is perceptually free from the effects of aliasing distortion when implemented at a sampling rate of F s = 88.2 kHz.A thorough evaluation of the proposed Lockhart model indicates that this configuration yields a signal quality equivalent to that of oversampling by a factor of eight (i.e., F s = 352.8kHz) at nearly a fourth of the computational expenses.For the case of the Serge wavefolder, the use of the antiderivative method produces an increase in signal quality equivalent to that of oversampling by a factor of two (i.e., F s = 88.2 kHz).
Furthermore, a recommended synthesis topology built around the Lockhart model consisting of four cascaded wavefolding stages, a saturator and a lowpass filter was presented.This topology was compared against a model of the Serge middle VCM built using six wavefolding stages.These structures illustrate the capabilities of wavefolding in a synthesis environment.However, it should be noted that the discussed topologies are not unique, as they can be modified according to the needs of the particular application.This effectively showcases the flexibility of VA models.

Supplementary Materials:
The following are available online at http://research.spa.aalto.fi/publications/papers/smc17-wavefolder: a real-time Max/MSP demo of the proposed Lockhart wavefolder topology implemented using Gen~and a MATLAB implementation of the Lambert-W function using Fritsch's iteration.

Figure 2 .Figure 3 .
Figure 2. Transfer function of the Lockhart wavefolder simulated using: (a) SPICE (Simulation Program with Integrated Circuit Emphasis); and (b) the proposed virtual analog (VA) model.Different colors indicate different values of R L .

Figure 4 .
Figure 4. Time-domain view of the proposed Lockhart wavefolder model plotted against its SPICE simulation for a 500-Hz sinusoidal input (peak amplitude 1 V) with load resistance: (a) R L = 10 kΩ; and (b) R L = 50 kΩ.

Figure 5 .
Figure 5. Schematic of a single folding cell in the middle section of the Serge Wave Multipliers (VCM).Figure adapted from [63].
Figure 5. Schematic of a single folding cell in the middle section of the Serge Wave Multipliers (VCM).Figure adapted from [63].

Figure 7 .
Figure 7. (a) Transfer function of a single wavefolding stage in the Serge middle VCM measured using SPICE and the proposed model; and (b) the absolute difference between these two curves.

Figure 8 .
Figure 8. Time-domain view of the Serge wavefolder model plotted against its SPICE simulation for a 500-Hz sinusoidal input with peak amplitude if 1 V.

FrequencyFigure 11 .
Figure 11.Spectrogram of the Lockhart wavefolder under 150-Hz sinusoidal input for values of R L between 1 and 50 kΩ.

FrequencyFigure 12 .
Figure 12.Spectogram for a linear sweep from 20 Hz to 5 kHz processed using: (a) the proposed Lockhart wavefolder model (R L = 50 kΩ); and (b) the proposed Serge wavefolder model.A sample rate F s = 1 MHz was used to simulate analog behavior.

FrequencyFigure 13 .
Figure 13.Spectrogram for a 1 V linear sweep from 20 Hz-5 kHz processed with the proposed Lockhart wavefolder (R L = 50 kΩ) model: (a) at audio rate; (b) using two times oversampling; (c) with antialiasing at audio rate; and (d) with antialiasing and oversampling by two.

FrequencyFigure 14 .
Figure 14.Spectrogram for a 1 V linear sweep from 20 Hz-5 kHz processed with the proposed Serge wavefolder model: (a) at audio rate; (b) using two times oversampling; (c) with antialiasing at audio rate; and (d) with antialiasing and oversampling by two.

Figure 15 .
Figure 15.Measured A-weighted noise-to-mask ratios (ANMRs) for a range of sinusoidal waveforms processed: (a) using the Lockhart wavefolder model (R L = 50 kΩ) under six different sampling rates; and (b) using the Serge wavefolder model under two different sampling rates, with and without the proposed antialiasing method.Values below the −10 dB threshold indicate lack of perceivable aliasing.

Figure 17 .
Figure 17.(a) Transfer function of the proposed Serge middle VCM; and (b) its output when driven by a 100-Hz sinewave for G S = 6 and zero dc offset.

FrequencyFigure 18 .
Figure 18.Spectrogram of: (a) a 150-Hz sinewave with amplitude 1 V processed by the proposed Serge middle VCM with varying gain G S from 0-6; and (b) a 200-Hz sinewave processed with varying gain G S from 0-3 and dc offset from 0-3 V.

Figure 20 .
Figure 20.(a) Transfer function of the proposed cascaded Lockhart wavefolder structure measured after the post-gain block; and (b) its output when driven by a 100-Hz sinusoidal input for G L = 10 and zero dc offset.

FrequencyFigure 21 .
Figure 21.Spectrogram of: (a) a 150-Hz sinewave peak amplitude 1 V processed by the proposed cascaded Lockhart topology with varying gain G L from 0 to15; and (b) a 200-Hz sinewave processed with varying gain G S from 0 to 0 and dc offset from 0 to 5 V.

Table 1 .
Component values for the Lockhart wavefolder circuit.

Table 3 .
Component information for the Serge wavefolder circuit shown in Figure5.

Table 4 .
Simulation parameters for a single folding stage in the middle Serge Wave Multiplier.

Table 5 .
Averaged computation times (in milliseconds) for the proposed Lockhart wavefolder model (R L = 50 kΩ) implemented in C for a 1-s 100-Hz sinusoidal input sampled at different oversampling (OS) rates and with different peak amplitude levels.