About a decade ago, I wrote two articles:
 Modulation Alternatives for the Software Engineer (November 2011)
 Isolated SigmaDelta Modulators, Rah Rah Rah! (April 2013)
Each of these are about deltasigma modulation, but they’re short and sweet, and not very indepth. And the 2013 article was really more about analogtodigital converters. So we’re going to revisit the subject, this time with a lot more technical depth — in fact, I’ve had to split this article into three parts — and a focus on the digitaltoanalog conversion side. The DAC case is interesting because all we’re doing is creating a series of low/high pulses at each clock tick, then putting them through a simple linear filter. There are very few opportunities for the introduction of unmodeled dynamics or sources of error, so we should expect very close results to any simulation.
We’ll be covering a number of aspects of firstorder and secondorder deltasigma modulators, with roughly the following outline:
 Part 1: Basics of modulation: kernels, block diagrams, linear analysis, choice of coefficients, and frequency spectra
 Part 2: Filtering tradeoffs, error analysis, and effects of low duty cycles
 Part 3: Implementation issues, including integer math and accumulator ranges
I’ll be using some theory, but you can learn a lot from empirical investigation; we’ll do the usual messing around in Python and see what we can find out.
Oh — and there’s one teensy little practical snag we’ll have to deal with. Please read all three parts, because I won’t bring it up until the discussion of implementation issues, and if you build a system based on the content of Part 1 or Part 2, you’re likely to run into this and it will ruin your day.
For now, though, let’s just keep it simple.
Author’s Note: I originally wrote this series of articles in January and February 2019, but put them on hold for various reasons. At the time, Parts 1 and 2 were essentially complete, and Part 3 was unfinished. I will try to publish Part 2 and Part 3 in a timely fashion, but it may take a while, so if you are reading this, please be patient.
Also a note about terminology: both “sigmadelta” and “deltasigma” are used in practice, and mean the same thing, but “deltasigma” seems to have the broader use. (See for example recent content from Texas Instruments and Analog Devices, and various books.)
Anyway: here we go!
FirstOrder DeltaSigma
Let’s start out with a Python implementation of the algorithm I mentioned in 2011, which has the following “kernel” to process an output \( y[k] \) from an input \( x[k] \) that is between 0 and 1:
accum1[k] = accum1[k1] + R*(x[k]  A*y[k1])
y[k] = 1 if accum1[k] >= R*A else 0
I’m going to use the numpy
and scipy
and matplotlib
packages.
%matplotlib inline import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import matplotlib.ticker as ticker import matplotlib import scipy.signal matplotlib.rcParams['mathtext.fontset']='cm'
def dsmod1(x, A, R=1, modulate=True): """ firstorder deltasigma modulation """ y = np.zeros_like(x) n = len(y) accum = 0 out = 0 M = 1 if modulate else 0 for i in range(n): accum += R*(x[i]  out*A) out_q = 1 if accum >= R*A else 0 out = (M*out_q +(1M)*accum) y[i] = out return y dsmod1.argformat = '' def samplewave(T=1,dt=0.002): t = np.arange(0,T+1e5,dt) x = 1np.abs(4*t1) x[4*t>2] = 0.95 x[t>0.75] = 0.02 return t,x def describe_dsmod(dsmod, args): sdtype = dsmod.__doc__.splitlines()[0].strip() try: sdargfmt = dsmod.argformat except: sdargfmt = '' return sdtype + sdargfmt.format(*args) def show_dsmod_samplewave(t,x,dsmod,args=(1,),tau=0.05, R=1, fig=None, return_handles=False, filter_dsmod=False): dt = t[1]t[0] if filter_dsmod: x1 = dsmod(x, *args,R=R, modulate=False) else: x1 = x y = dsmod(x,*args,R=R) a = dt/tau yfilt = scipy.signal.lfilter([a],[1,a1],y) xfilt = scipy.signal.lfilter([a],[1,a1],x1) if fig is None: fig = plt.figure(figsize=(10,7), facecolor='white') gs = gridspec.GridSpec(3, 1,height_ratios=[2,0.4,1]) axes = [fig.add_subplot(gs[k]) for k in range(3)] ax0 = axes[0] ax0.plot(t,x,label='x(t)') ax0.plot(t,yfilt,label='LPF(y(t))') ax0.plot(t,xfilt,label='LPF(x(t))') ax0.legend(labelspacing=0,loc='upper right',fontsize=10) ax0.set_ylabel('x, LPF(x), LPF(y)') ax0.set_xticklabels([]) ax1 = axes[1] ax1.step(t,y,linewidth=0.5) ax1.set_ylim([0.1,1.1]) ax1.set_yticks([0,1]) ax1.set_ylabel('y') ax1.set_xticklabels([]) ax2 = axes[2] filt_error = xfiltyfilt ylim = (min(filt_error),max(filt_error)) # plot margin s = 0.08 ylim = (ylim[0]*(1+s)ylim[1]*s, ylim[1]*(1+s)ylim[0]*s) ax2.plot(t,filt_error) ax2.set_ylim(ylim) for s in [1,1]: # approximate error bounds; we'll show why later in the article ax2.plot([t[0],t[1]],[s*a/2,s*a/2],color='#008000',linewidth=0.8,dashes=[8,8]) ax2.set_ylabel('error') ax2.get_yaxis().set_major_locator(ticker.MaxNLocator(12,steps=[1,2,5,10])) ax2.set_xlabel('time (s)') ax2.grid(True,linestyle=':') for ax in axes: ax.set_xlim(min(t),max(t)) fig.subplots_adjust(hspace=0.1) fig.suptitle('Deltasigma modulation (%s)\nwith RC filter $\\tau=%.2f s = %.1f\\Delta t, R=%.2f$' % (describe_dsmod(dsmod,args), tau, tau/dt, R)) if return_handles: return fig, (ax0,ax1,ax2) t,x = samplewave() show_dsmod_samplewave(t,x,dsmod1,filter_dsmod=True)
Now, there is an amazing amount of information we can see in the graph above, from just a handful of lines of Python. What are we looking at?
The top subplot contains \( x(t) \), a sample waveform that starts with a triangle wave and ends with a square wave. We feed this into our 1storder deltasigma modulator and we get out a pulse train \( y(t) \) in the middle subplot, which sort of looks like a PWM waveform except that it isn’t constant frequency. We have an accumulator, and whenever it gets to some threshold \( R\times A \), we output a 1, and next time we subtract \( R\times A \) from the accumulator; otherwise we output a zero.
The DC gain from input to output is \( 1/A \). If you have floatingpoint data limited between 0 and 1, use \( A=1 \). This will produce unity gain at low frequencies; since we’re integrating \( x(t)Ay(t) \), if \( A=1 \) then the lowfrequency components of \( x(t) \) and \( y(t) \) balance each other out. If you have fixedpoint integer data limited between 0 and some value \( A \), use that \( A \) instead.
The fun comes when we filter the input \( x(t) \) and the output \( y(t) \) with a firstorder lowpass filter. The results are shown in the top subplot, and they follow each other very well, because at low frequencies \( x(t) \) and \( y(t) \) are essentially identical. The error between the filtered waveforms is shown in the bottom subplot, and it’s essentially limited to a small band that turns out to be limited to approximately \( \pm\Delta t/2\tau \) for smoothlychanging inputs, as we’ll see later. (In this case, \( \Delta t = 0.002 \)s and \( \tau=0.05 \)s, so the theoretical value of error band amplitude is \( \pm0.02 \), which matches our example waveform fairly well.)
This is a very simple DAC; we implement our deltasigma modulator digitally, and take our pulse waveform and feed it through an RC filter.
Linear Analysis of 1storder DeltaSigma Modulation
The modulator itself is usually analyzed with a linear approximation for the threshold comparator, namely that we replace it with an additive noise waveform \( Q(t) \), or \( Q[k] \) for discretetime analysis, that is uncorrelated with anything else. (There are only two ways we can get away with this. One is if the integrator output is zooming around so much that the difference between the integrator output and the modulator output looks enough like noise. This happens if the input frequency spectrum is “interesting” enough, or in certain cases of “uninteresting” input, like DC inputs where the DC level is not near 0 or 1. The second possibility is if we add an intentional “dither” noise signal to the threshold comparator input, improving the validity of this noise approximation, at the expense of additional output noise.)
Here \( \Delta t \) is the timestep for the discretetime system; \( T_0 = \Delta t / R \) is a tunable integrator parameter, but it’s pretty standard to use \( T_0 = \Delta t \) so that \( R = \Delta t/T_0 = 1. \)
The transfer function from input \( x(t) \) to output \( y(t) \) is \( H(z) = R/(1z^{1})/(1+Rz^{1}/(1z^{1})) = \dfrac{R}{1+Rz^{1}  z^{1}} \) which has a pole at \( z=1R \), and is therefore stable for \( 0 < R < 2 \), with marginal stability when \( 0 < R \ll 1. \) The noise transfer function from quantization noise \( Q(t) \) to output \( y(t) \) is \( N(z) = 1/(1+Rz^{1}/(1z^{1})) = \dfrac{1z^{1}}{1+(R1)z^{1}} \).
Let’s graph these for a few values of \( R \):
fig = plt.figure(figsize=(8,6), facecolor='white') ax = fig.add_subplot(1,1,1) f = np.arange(0,0.5,0.001) # frequency, normalized to 1/delta T for R in [0.1,0.5,1.0,1.5]: z = np.exp(2j*np.pi*f) zinv = 1/z H = R/(1+(R1)*zinv) N = (1zinv)/(1+(R1)*zinv) line = ax.plot(f,np.abs(H), label='$H(z), R=%.1f$' % R) linecolor = line[0].get_color() ax.plot(f,np.abs(N),color=linecolor,linestyle='', label='$N(z), R=%.1f$' % R) x0 = 0.5 ax.plot([0,x0],[0,x0*(2*np.pi/R)],'',color=linecolor,linewidth=0.4) ax.legend(labelspacing=0, fontsize=12, loc='upper left') ax.set_xlabel('$f\Delta t$', fontsize=16) ax.set_xlim(0,0.5) ax.set_ylabel('$H(z), N(z)$', fontsize=16) ax.set_ylim(0,4.1);
This shows that \( H(z) \) is a filter with unity gain at low frequencies — for \( R<1 \) it is a lowpass filter; for \( R>1 \) the gain actually increases. On the other hand, its complement \( N(z) \) suppresses frequency content at low frequencies, where \( N(z) \approx 2\pi f\Delta t/R \). This is the intended effect; we want the remaining noise energy to be in the higher frequency ranges so that a lowpassfilter on the output will get rid of most of it.
(I should note in passing that I dislike using ztransform analysis, mostly because of my lack of familiarity. It’s very useful, but I find it a lot less intuitive than Laplace (s) domain analysis, so sometimes I just pretend that a system is continuoustime and use sdomain analysis instead. This works okay as long as the frequency content is small once it gets to near Nyquist frequency \( f_s/2 \), otherwise you need ztransforms to stay accurate.)
With linear analysis, it looks like the value of \( R \) could be increased to help attenuate quantization noise further. Unfortunately, this is not actually the case; the assumption of uncorrelated noise breaks down. We can try looking at the effect of changing the integrator gain \( R \) in a test simulation of a deltasigma modulator, but it doesn’t change the behavior at all. This seems really weird, until you look more closely at the kernel equations:
accum1[k] = accum1[k1] + R*(x[k]  A*y[k1])
y[k] = 1 if accum1[k] >= R*A else 0
We can make the substitution accum1[k] = R*acc_n[k]
, at which point we have
acc_n[k] = acc_n[k1] + x[k]  A*y[k1]
y[k] = 1 if acc_n[k] >= A else 0
This is the same behavior, independent of \( R \), so all that a change of \( R \) is doing is to change the scaling of the integrator value relative to input \( x \) and output \( y \).
Below are two other graphs with \( R=0.01 \) and \( R=100.0 \).
t,x = samplewave() show_dsmod_samplewave(t,x,dsmod1,R=0.01)
t,x = samplewave() show_dsmod_samplewave(t,x,dsmod1,R=100.0)
See? Changing \( R \) doesn’t make any difference. So we’ll just use \( R=1 \) from here onward.
Frequency Spectrum of 1storder DeltaSigma Modulation
Let’s get a little more sophisticated and measure the frequency spectra of the modulated output from a couple of input waveforms \( x(a_0,a_1,t) = a_0 + a_1 \sin \omega_1 t, \) which is just a sinewave plus a DC offset.
def enumerate_reverse(somelist): n = len(somelist) while n > 0: n = 1 yield n, somelist[n] def f1_noise(t,a0,a1): r = np.random.RandomState(seed=123) return (r.randint(2,size=t.shape)*21)*a1 + a0 def freq_spectrum_sd(dsmod, args=(), N=65536,T=1.0,Tsettle=0.1, f1=16,a1=0.01,a0list=(0.5,0.4,0.3,0.2,0.1,0.05,0.03,0.02,0), dBmin = 80, window=None, envelope_kmax=0, envelope_min=1, zreverse=False, return_handles=False, fig=None, show_time_waveforms=True): """ f1: frequency of signal to modulate, *OR* a callable: f1(t,a0,a1) Tsettle: settling time before displaying or computing FFT """ dt = T/N Nsettle = int(np.floor(Tsettle/dt)) t = np.arange(Nsettle, N)*dt f = np.arange(N)/T if window is None: W = np.ones(N) else: W = window(N) if fig is None: fig = plt.figure(figsize=(10,12)) nwaveforms = 3 if show_time_waveforms else 2 if show_time_waveforms: ax0 = fig.add_subplot(nwaveforms,1,1) ax0.set_xlabel('time (s)') ax0.grid(True) ax0.set_xlim(0,T) ix1 = 2 ix2 = 3 else: ax0 = None ix1 = 1 ix2 = 2 ax1 = fig.add_subplot(nwaveforms,1,ix1) ax2 = fig.add_subplot(nwaveforms,1,ix2) fig.subplots_adjust(hspace=0.25) for ax,fmax in [(ax1,32768),(ax2,2000)]: ax.set_xlabel('freq (Hz)') ax.set_ylabel('amplitude (dB)') ax.set_ylim([dBmin,0]) ax.set_xlim([0.01*fmax,fmax]) ax.set_yticks(np.arange(dBmin,0,5), minor=True) ax.grid(True) ax.grid(which='minor', linestyle=':') colors = ['blue','green','red','turquoise','magenta','orange','black','skyblue','brown'] e = enumerate_reverse if zreverse else enumerate zero_margin = 0.0002 for i,a0 in e(a0list): a0 = max(a0,a1+zero_margin) if callable(f1): x = f1(t,a0,a1) else: x = a0+a1*np.sin(2*np.pi*f1*t) y = dsmod(x, *args) Y = np.fft.fft(y[Nsettle:]*W)/N Y[1:] *= 2 # all except DC have half their energy in positive and negative frequencies dBY = 20*np.log10(np.abs(Y)+1e12) dBY0 = dBY.copy() # find the maximum over 2^envelope_kmax points dBY[:envelope_min] = 1000 for k in range(envelope_kmax): nshift = 1<<k dBY[nshift:] = np.maximum(dBY[nshift:], dBY[:nshift]) dBY[:envelope_min] = dBY0[:envelope_min] color=colors[i%len(colors)] if show_time_waveforms: ax0.plot(t[Nsettle:],x[Nsettle:],color=color) ax1.plot(f,dBY,color=color) ax2.plot(f,dBY,color=color) fig.suptitle('Frequency spectrum of unfiltered modulated signals,\n' +describe_dsmod(dsmod,args), y=0.92) if return_handles: return fig, (ax0,ax1,ax2) freq_spectrum_sd(dsmod1,[1.0])
Again, a world of information here in just a simple example. The top plot is the timeseries input waveforms. The bottom two plots are the FFT magnitude in decibels, at different frequency scales, and they contain three components:
 the DC component at \( f=0 \)
 the sinewave component: 0.02 amplitude = 40dB at \( f=16 \) Hz
 modulation noise
The interesting thing here, of course, is the modulation noise. We can’t get rid of it, but the frequency spectrum depends very strongly on the DC offset. In all cases, there’s very little frequency content at 1kHz and below, with the lowest amplitudes nearest to DC. Beyond that, there are different ranges where the quantization noise kicks in, and it depends upon the DC offset. When the DC offset is near 0.5, we get quantization noise near \( f_s/2 \), which makes sense: here the deltasigma output is chattering back and forth between 0 and 1 on average every other sample. Lower DC offsets have the chattering at lower frequencies since, for example, at a DC offset of 0.2, on average out of every 5 samples, one of them will be a 1 and the rest will be a zero, so we could expect \( f_s/5 \approx 13107 \) Hz in this case. (The same is true for DC offsets greater than 0.5; we haven’t shown them here, but you can use symmetry arguments to show that an input of \( 1x(t) \) will have identical output to that for \( x(t) \) but with the 0’s and 1’s swapped.)
Note that this is the spectrum of the modulated signals. We haven’t done any filtering, we’re just looking at what a bunch of cleverlygenerated pulses is like in the frequency domain.
The frequency spectra (in the middle graph above) have “ears” because the inputs \( x(t) \) are not constant, but are slowlyvarying with respect to the deltasigma output frequency. More specifically, inputs \( x(t) \) are sinusoidal and spend more time at the extremes; we should expect the frequency shape to look similar to a histogram of a sine wave, plotted in log scale:
N = 65536 t = np.arange(N)*1.0/N a1 = 0.01 w1 = 2*np.pi*16 x = 0.4+a1*np.sin(w1*t) hcount,hx = np.histogram(x,100) hx_centers = (hx[:1]+hx[1:])/2.0 fig,ax = plt.subplots(figsize=(8,6), facecolor='white') ax.plot(hx_centers,20*np.log10(hcount)) ax.get_yaxis().set_major_locator(ticker.MaxNLocator(12,steps=[1,2,5,10])) ax.set_ylabel('Relative frequency, dB') ax.grid(True);
In this case, the “ears” are much larger — about a 20 dB discernible difference between the extremes and the center, whereas in the deltasigmamodulated spectra, the extremes are about 45 dB higher. But it’s the same general effect.
The individual frequency spectra are easier to see if we plot the spectral envelope, that is, the maximum value over some frequency band, so that the plots don’t obscure each other so easily. The envelope_kmax
and envelope_min
arguments can help; here’s the same spectra but with an envelope over \( 2^6 = 64 \) successive frequency points (envelope_kmax=6
) starting at a minimum of 50Hz:
freq_spectrum_sd(dsmod1,[1.0],envelope_kmax=6, envelope_min=50, dBmin=100)
We can see the spectral distribution of quantization noise fairly easily. In all cases the noise content is lower at low frequencies, and doesn’t really reach 50dB until somewhere between 7006000Hz depending on the DC offset. We don’t care about the noise content at high frequencies since the intent is to attenuate it with a lowpass filter.
Now, the interesting thing is that when you use a firstorder lowpass filter — which is just a plain old RC filter in analog — the error between the filtered output and the filtered input (\( =x^F(t)  y^F(t) \) where I’m using this superscript F to denote a lowpass filter) has a bounded amplitude. The easiest way to think about it, is that at each timestep \( \Delta t \), the filtered output \( y^F(t) \) can take a trajectory upwards with a 1 output, or downwards with a 0 output; the filtered input \( x^F(t) \) is somewhere in between, and the difference between these possible outputs is bounded. We’ll do the math in a second, but here’s what I’m talking about:
a0 = 0.4 di = 0.01 dt_tau = 0.05 # larger than the example above, to show curvature Ni = 20 ilist = np.arange(0,Ni,di) fig = plt.figure(figsize=(9,6)) ax = fig.add_subplot(1,1,1) ax.plot(0,a0,'.') ai = [a0,a0] for ii in range(6): iifuture = ilist[ilist>=ii] iistep = iifuture[iifuture<=ii+1] f0 = lambda i,a: a*np.exp((iii)*dt_tau) f1 = lambda i,a: 1(1a)*np.exp((iii)*dt_tau) for f,c,w in [(f0,'red',0),(f1,'green',1)]: a = ai[1w] ax.plot(iifuture,f(iifuture,a),color=c,linestyle=':') ax.plot(iistep,f(iistep,a),color=c) ax.plot(ii+1,f(ii+1,a),'.',color=c) ai = [f0(ii+1,ai[1]), f1(ii+1,ai[0])] ax.set_xlim(Ni*0.01,Ni) ax.get_xaxis().set_major_locator(ticker.MaxNLocator(12,steps=[1,2,5,10])) ax.set_xlabel('# steps') ax.set_ylabel('Filtered output $y^F$') ax.set_title('Potential trajectories of filtered deltasigma modulation, firstorder lowpass filter');
At each step, the possible trajectories form a sort of a herringbone pattern, with 1 outputs forming the upward curves \( y^F = 1(1y^F_k)e^{(tt_k)/\tau} \) and 0 outputs forming the downward curves \( y^F=y^F_ke^{(tt_k)/\tau} \). If we let \( \rho = e^{\Delta t/\tau} \), then the upwards step change is \( \Delta y^F_+ = 1\rho(1y^F_k) \) and the downwards step change is \( \Delta y^F_ = \rho y^F_k \), and the span between them is \( \Delta y^F_+  \Delta y^F_ = 1\rho(1y^F_k)  \rho y^F_k = (1\rho). \) This is a constant, independent of the filtered output \( y^F \), and it represents the interval in which the filtered input \( x^F \) lies, so the maximum error magnitude is just this interval length \( 1\rho = 1e^{\Delta t/\tau} \approx \Delta t/\tau. \) It turns out that for slowlychanging inputs, the maximum error is about half of this, or \( \approx \Delta t/2\tau. \)
If we look at the first graph I showed in the article, we have \( \Delta t = 0.002 \) and \( \tau=0.05 \), so \( \Delta t/\tau=0.04 \) and half of this is 0.02, which is the error bound amplitude shown in the graph.
Secondorder DeltaSigma Modulation
Things get a little more interesting when we go to secondorder deltasigma modulation. Let’s look at the same kind of error graphs and frequency spectra as in the firstorder deltasigma modulation case.
(We’re going to make one change though: to compare filtered input vs. filtered \( \Delta\Sigma \) modulation, we’ll feed the input into an “unmodulated” version of the modulator, which captures the transfer function \( H(z) \)’s lowpass filtering effect without the quantization.)
Here’s the kernel:
accum1[k] = accum1[k1] + R*(x[k]  A*y[k1])
accum2[k] = accum2[k1] + R*(accum1[k]  B*y[k1])
y[k] = 1 if accum2[k] >= R*B else 0
We have two accumulators, and two independent feedback gains \( A \) and \( B \).
def dsmod2(x, A, B, ffA=0, ffB=0, ffC=0, R=1, modulate=True): """ 2ndorder deltasigma modulation ffA, ffB, ffC are feedforward terms not normally used """ y = np.zeros_like(x) n = len(y) accum1 = 0 accum2 = 0 y_out = 0 M = 1 if modulate else 0 for i in range(n): accum1 += R*(x[i]  A*y_out) accum2 += R*(accum1 + x[i]*ffA  B*y_out) q_in = accum2 + x[i]*ffB + accum1*ffC y_out = (M*(1 if q_in >= R*B else 0) +(1M)*q_in) y[i] = y_out return y dsmod2.argformat = ', $B={1:.3f}$' t,x = samplewave() fig,axes=show_dsmod_samplewave(t,x,dsmod2,args=(1,1.0), return_handles=True,filter_dsmod=True) import matplotlib.patches as mpatches for ax in axes: yl = ax.get_ylim() ax.add_patch(mpatches.Rectangle((0.49,yl[0]),0.27,yl[1]yl[0], fill=True, color='orange', alpha=0.12, figure=fig)) axes[2].annotate('', xy=(0.6, 1.3), xycoords='axes fraction', xytext=(22, 35), textcoords='offset points', arrowprops=dict(arrowstyle=">", connectionstyle="arc3", shrinkA=0, ), ) axes[2].annotate('Note "galloping"', fontsize=12, xy=(0.6, 0.03), xycoords='data', xytext=(20, 12), textcoords='offset points', bbox=dict(boxstyle="round", fc="w"), arrowprops=dict(arrowstyle=">", connectionstyle="arc3", shrinkA=0, relpos=(0.1,0.5)), );
freq_spectrum_sd(dsmod2,(1.0,1.0),envelope_kmax=6, envelope_min=50, dBmin=100)
This looks mostly the same as the 1storder graphs, with a few notable exceptions:
 The noiseshaping at high frequencies (2kHz and higher) is a little better. We can see this in two ways:
 The peak “ears” are 2–5 dB lower with secondorder modulation
 Secondary ears are essentially gone in secondorder modulation — the remainder of the spectrum is smoother, spreaded more equally as noise.
 Noiseshaping at midrange frequences (200Hz  2kHz) is better in most cases, with the exceptions at low DC offset.
 For DC offsets of 0.1 and above, with secondorder modulation, the frequency content dropped below a curve leading roughly from 88dB at 500Hz to 66dB at 2kHz.
 With the 1storder deltasigma modulation, frequency content of quantization noise at 500Hz and below was in the 80dB to 57dB range; here for the 2ndorder deltasigma modulator it’s below 88dB except for the very lowest DC offsets (85dB at \( a_0=0.03 \), 72dB at \( a_0=0.02 \), 57dB at \( a_0=0.0102 \)). The only example that got worse below 500Hz was \( a_0 = 0.02 \), the light blue curve, which increased in the 250Hz – 600Hz range.
 At low frequencies (below 1000Hz) due to DC offsets close to 0 or 1, the error after filtering is worse with 2ndorder deltasigma modulation than with 1storder deltasigma modulation.
 This is visible between 500Hz – 1250Hz at \( a_0=0.03 \) (black curve), 250 – 600Hz at \( a_0 = 0.02 \) (light blue), and below 200Hz at \( a=0=0.0102 \) (brown).
 In part, this is caused by the presence of subharmonic idle tones; note the “galloping” effect in the output and error timeseries plots, where the durations between pulses are of alternating lengths (long, short, long, short). We’ll talk about this a little bit later, but one way of thinking about it is that secondorder deltasigma modulation is more “clever” than firstorder deltasigma modulation; the cleverness works well as long as the duty cycle doesn’t get really close to 0 or 1, where the integrators take more time to make changes in the output.
SecondOrder DeltaSigma Modulation Coefficients
One obvious question to ask is where do these A, B coefficients come from, for the 2ndorder deltasigma modulator. The A coefficient is the same as in the firstorder modulator (1 for floatingpoint values, maximum input for fixedpoint values); it determines the lowfrequency scaling factor. The B coefficient can be varied as a compromise between stability and quantization noise attenuation, but the A=1, B=1 values seem to be about the best in general, with larger values of B to cover the “galloping” cases where the inputs are close to the limits — again, we’ll talk about this later.
If you’re venturing into 3rdorder or higher deltasigma modulation, you need to know what you’re doing and get comfortable with ztransforms in order to ensure stability — but 1storder and 2ndorder deltasigma are fairly simple to get something workable.
If you really want to get picky, we can analyze the signal and noise transfer functions of a 2ndorder deltasigma modulator for stability. (If not, just skip ahead past the math.) Remember that \( R=\Delta T/T_0 = 1 \) again, and in this case:
$$Y(z) = H(z)X(z) = R/(1z^{1})\left(R/(1z^{1})(X(z)  z^{1}H(z)X(z))  Bz^{1}H(z)X(z)\right)$$
$$H(z) = R/(1z^{1})\left(R/(1z^{1})(1  z^{1}H(z))  Bz^{1}H(z)\right)$$
$$H(z) + R^2z^{1}/(1z^{1})^2H(z) + BRz^{1}/(1z^{1})H(z) = R^2/(1z^{1})^2$$
$$\begin{align}H(z) &= \frac{R^2/(1z^{1})^2}{1 + R^2z^{1}/(1z^{1})^2 + BRz^{1}/(1z^{1})} \cr &= \frac{R^2}{(1z^{1})^2 + R^2z^{1} + BRz^{1}(1z^{1})} \cr &= \frac{R^2}{12z^{1}+z^{2} + R^2z^{1} + BRz^{1}  BRz^{2}} \cr &= \frac{R^2}{1 +(2+R^2+BR)z^{1}+(1BR)z^{2}} \end{align}$$
$$\begin{align}N(z) &= \frac{1}{1 + R^2z^{1}/(1z^{1})^2 + BRz^{1}/(1z^{1})} \cr &= \frac{(1z^{1})^2}{1 +(2+R^2+BR)z^{1}+(1BR)z^{2}} \end{align}$$
which both have poles at the zeros of \( 1 +(2+R^2+BR)z^{1}+(1BR)z^{2} \), which occurs when \( z^2 + (2+R^2+BR)z +(1BR) = 0 \) or in other words at
$$\begin{align}z &=\frac{2R^2BR \pm \sqrt{(2R^2BR)^2  4(1BR)}}{2} \cr &=1(R^2+BR)/2 \pm \sqrt{(1(R^2+BR)/2)^2  1+BR} \cr &=1+R/2\left((R+B) \pm \sqrt{(R+B)^24}\right). \cr \end{align}$$
This is stable if the poles occur at \( z < 1 \).
If \( 0 \le (R+B) < 2 \) then we have a pair of complex poles with \( z^2 = (1R^2/2BR/2)^2 + R^2/4(4(R+B)^2) = 1BR \) which is between 0 and 1 as long as \( B \) and \( R \) are both positive.
If \( (R+B) \ge 2 \) then we have real poles at \( z = 1+R/2((R+B) \pm \sqrt{(R+B)^2  4}) \). The mess in parentheses is always negative, regardless of the plus/minus sign.
My skills at Grungy Algebra are fading at the moment, so let’s just shown what happens with \( R=1 \). In this case we have \( B>1 \) which produces real poles at \( z = 1(1+B)/2 \pm \sqrt{(1+B)^24}/2 = (1B)/2 \pm \sqrt{(1+B)^24}/2 \). Both poles are within the unit circle as long as \( B<1.5 \); at \( B=1.5 \), the poles are at \( z=1.0, +0.5 \).
So for the \( R=1 \) case, the system should be stable as long as \( 0 < B < 1.5 \).
At least, that’s the theory, and here are plots of the input and noise transfer functions:
fig = plt.figure(figsize=(9,7)) ax = fig.add_subplot(1,1,1) f = np.arange(0,0.5,0.001) # frequency, normalized to 1/delta T R = 1 for B in [0.5,0.6,0.8,1,1.2]: z = np.exp(2j*np.pi*f) zinv = 1/z denom = 1+(2+R*R+B*R)*zinv+(1B*R)*zinv*zinv H = R*R/denom N = (1zinv)**2/denom line = ax.plot(f,np.abs(H), label='$H(z), B=%.1f$' % B) ax.plot(f,np.abs(N),color=line[0].get_color(),linestyle='', label='$N(z), B=%.1f$' % B) ax.legend(labelspacing=0, fontsize=12, loc='upper left') ax.set_xlabel('$f\Delta t$', fontsize=16) ax.set_ylabel('$H(z), N(z)$', fontsize=16) ax.set_ylim(0,4) ax.set_xlim(0,0.5) ax.set_title('Input transfer function $H(z)$, noise transfer function $N(z)$', fontsize=16);
The \( N(z) \) curves at low values of \( f\Delta t \) are almost right on top of each other; they approach the envelope \( \leftN(z)\right = (\omega\Delta t)^2 = (2\pi f\Delta t)^2 \):
fig = plt.figure(figsize=(9,7)) ax = fig.add_subplot(1,1,1) f = np.arange(0,0.11,0.0005) # frequency, normalized to 1/delta T R = 1 for B in [0.5,0.6,0.8,1,1.25]: z = np.exp(2j*np.pi*f) zinv = 1/z denom = 1+(2+R*R+B*R)*zinv+(1B*R)*zinv*zinv N = (1zinv)**2/denom ax.plot(f,np.abs(N), label='$N(z), B=%.1f$' % B,linewidth=0.5, linestyle='' if B == 2 else '') ax.plot(f,(2*np.pi*f)**2,color='black',linestyle='',linewidth=1, label='limit: $(2\pi f\Delta t)^2$') ax.legend(labelspacing=0, fontsize=10, loc='upper left') ax.set_xlabel('$f\Delta t$', fontsize=13) ax.set_ylabel('$H(z), N(z)$', fontsize=13) ax.set_xlim(0,0.1) ax.set_ylim(0,0.5);
In practice, things are a little different than a linearized model might suggest. Low values of \( B \) can cause large transients.
In the graphs below we look at several values of \( B \), and for calculating the output error we include the transfer function of the linearized (nonmodulating) system in the \( LPF(x(t)) \) path to equalize the effects of time delay.
t,x = samplewave() for B in [0.4, 0.6, 0.8, 1.0, 1.1, 1.2, 1.3, 1.4, 1.49]: show_dsmod_samplewave(t,x,dsmod2,args=(1,B),filter_dsmod=True)
This is an informal demonstration that \( B \) in the 1.2  1.4 range seems to be the best at minimizing error and maintaining stability. (Near 1.5, the system is marginally stable and it takes a very long time for oscillations to die out.) If the modulator input is a DC value close to the limits, then lower values of \( B \) allow a “galloping” effect in the output, producing subharmonic idle tones.
Wrapup
OK, so what have we learned today?

The “kernel” (defining equation) of a first and secondorder modulator is very simple:

firstorder:
accum1[k] = accum1[k1] + R*(x[k]  A*y[k1]) y[k] = 1 if accum1[k] >= A*R else 0

secondorder:
accum1[k] = accum1[k1] + R*(x[k]  A*y[k1]) accum2[k] = accum2[k1] + R*(accum1[k]  B*y[k1]) y[k] = 1 if accum2[k] >= B*R else 0


The corresponding block diagrams are shown below, and can be linearized by replacing the quantization (comparator) block with an additive source of quantization noise. (In the diagrams, \( R = \Delta t/T_0 \), and I have chosen \( A=1 \); the gain \( A \) is located in the negative input path of the first summing junction.)
 firstorder:
 secondorder:

The value of \( A \) controls the lowfrequency gain \( 1/A \) from input to output, and should usually be 1 in floatingpoint systems.

The value of \( R \) controls the integration time constant, but at least for firstorder systems, the value of \( R \) affects only the scaling of the integrator, so just use \( R=1 \).

The error band of a firstorder modulator with timestep \( \Delta t \), when passed through a singlepole RC filter with time constant \( \tau = RC \), is \( \pm\Delta t/2\tau \).

Linear analysis assumes that the output quantizer behaves like a source of additive uncorrelated noise, in which we can determine the inputoutput transfer function \( H(z) \) and the quantization noise transfer function \( N(z) \).

We can look at the frequency spectrum of modulated signals, and as long as the inputs aren’t too close to 0 or \( A \), then the noise gets suppressed at lower frequencies and pushed out to higher frequencies; the noise transfer function \( N(z) \) behaves roughly as \( \leftN(\omega)\right \approx \omega \Delta t \) for firstorder modulation and \( \leftN(\omega)\right \approx (\omega \Delta t)^2 \) for secondorder modulation.

For secondorder systems, \( B \) should be chosen in the 1.0  1.4 range, with higher values used to reduce a “galloping” behavior.
Next time we’ll look at some of the tradeoffs of output filters for deltasigma modulation.
© 2019, 2023 Jason M. Sachs, all rights reserved.