# How would sigma-delta modulation work in software?

I've recently come across this YouTube video that shows a PCM stream being played back via the 1-bit internal "PC Speaker", and after hours of googling I've come to think this is done via Sigma-Delta modulation.

So far, what I gather from the process is that it takes the original digital signal, oversamples and approximates it via PWM while applying a low-pass filter to the 1-bit output. I'm not quite sure if this is right.

Not having a background in DSP, I'm wondering:
How could this process be implemented in software?

1. Upsample the signal to desired frequency. You can use linear interpolation (poor quality but simple and fast) or sinc interpolation (good quality but slow). For example of library dedicated to good quality resampling, see libsoxr: http://sourceforge.net/p/soxr/wiki/Home/

2. Perform dithering with noise shaping. You need to filter the dithering signal to minimize noise within 0-20000 Hz range while increasing it in ultrasonic range. For details on noise shaping, see this excellent document: http://www.beis.de/Elektronik/DeltaSigma/DeltaSigma.html especially Figure 9 which shows adding noise in feedback loop.

That's how delta-sigma DACs work and this can certainly be implemented in software.

PWM would require upsampling 65536× for 16-bit audio so it isn't used. PDM (Pulse Density Modulation) is used instead and is essentially the same as delta-sigma modulation. And it requires only 64× upsampling for Hi-Fi audio.

• Is the dithering step before modulating the signal with PDM? Nov 5 '13 at 15:05
• Dither noise is added in PDM's feedback loop, so it's not before or after PDM, it's part of delta-sigma modulator.However, if noise shaping isn't needed (because you are converting, say, 32bit float to 16bit int) dither noise can be added before bit depth reduction, but it's a different story. Nov 5 '13 at 15:22

Here is a detailed implementation that may provide further insight. This is a 2nd order Sigma Delta DAC that I have implemented in Python. The simplified but completely functional code is below (finer details regarding bit and cycle accurate operation have been omitted for clarity).

Below is the functional block diagram. Noise shaping is specifically provided because of the feedback structure. (See my answer here for more details on that: How to plot noise shaped spectrum of First order Incremental Sigma Delta ADC's output? ). This is no different than any control loop which tracks a reference (low pass filter to reference input which is the desired level) but allows high frequency of the controlled device to pass through (high pass filter to the quantization noise).

Example Python Code

def dsigma(input_values, width):
# create output by iterating through input_values
for next_sample in input_values:

# compute next state (clock update)
sum1d = sum1
sum2d = sum2

# asynchronous operations
out = -1 if sum2d < 0 else 1
fb = -2**(width-1) if sum2d < 0 else 2**(width-1)-1
fbx2 = 2*fb
delta1 = next_sample - fb
delta2 = sum1d - fbx2

sum1 = sum1d + delta1
sum2 = sum2d + delta2
yield out


This is a python generator, so you could use list or in a loop to generate each value. For example below shows a simple case of 5 samples for a 24 bit input. (You would need to generate much longer sample set to see the Sigma Delta PWM operation in action):

out = dsigma([25,25,25,25,25], 24)
print(list(out))


The output spectrum from this simple example demonstrates the noise shaping which results in much higher effective precision from the two level output once filtered (eliminating all the high frequency noise and providing a precision output level representing the average of the 2 level PWM signal.)

• somewhere i have some cheesy MATLAB code for a 2nd-order sigma delta. no upsampling, but it's amazing how much of the music you can hear with 44.1 kHz 1-bit audio. May 2 '20 at 1:12
• 12 dB/octave so you get 2 bits for every doubling of the sampling rate—- so I guess you were getting the equivalent of 4 bits at 11 KHz sampling May 2 '20 at 1:15

hay i found my sigma delta MATLAB program. dunno if it's any good.

%
%
%
%
%
%
%   simulated 1 bit sigma-delta converter:
%
%
%            x(n)-y(n-1)    w(n)                v(n)                 ( mean(y^2) = A^2 )
%
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[Quantizer]----.---> y = +/- A = quantized value
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%   "linearized" model:
%                                                          .---- q = quantization noise  ( mean(q) = 0 )
%                                                          |
%                                                          |
%            x - y/z        w                   v          |         ( mean(y^2) = G^2*mean(v^2) + mean(q^2) )
%                                                          v
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[G]--->(+)-----.---> y = G*v + q
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%
%
%
%           W = 1/(z-1)*(X - Y/z)
%
%
%           V = 1/(z-1)*(W - fbg*Y/z)
%
%             = (X - Y/z - fbg*Y*(z-1)/z)/(z-1)^2
%
%             = (X*z - Y*(1+fbg*(z-1))) / (z*(z-1)^2)
%
%
%           Y = G*V + Q = G*(X*z - Y*(1+fbg*(z-1)))/(z*(z-1)^2) + Q
%
%             = G*X/(z-1)^2 - G*Y*(1+fbg*(z-1))/(z*(z-1)^2) + Q
%
%
%           Y + G*Y*(1-fbg + fbg*z)/(z*(z-1)^2) = G*X/(z-1)^2 + Q
%
%
%           Y = (G*X/(z-1)^2 + Q)/(1 + G*(1-fbg + fbg*z)/(z*(z-1)^2))
%
%             = (G*X/(z-1)^2 + Q)*(z*(z-1)^2)/((z*(z-1)^2) + G*(1-fbg + fbg*z))
%
%             = z*(G*X + Q*(z-1)^2)/(z^3 - 2*z^2 + (G*fbg+1)*z + G*(1-fbg))
%
%             = z*(G*X + Q*(z-1)^2)/(z*(z-1)^2 + G*fbg*z + G*(1-fbga))
%
%
%    as z -> 1  (DC)
%
%           Y  ->  z*X/(fbg*z + (1-fbg)) =  X/(fbg + (1-fbg)/z)  -->  X
%
%

if ~exist('mean_vv', 'var')
linearized_model = 0                % run this with 0 the first time to define G and mean(q^2)
end

if ~exist('A', 'var')
A = 1.0                             % comparator output magnitude
end

if ~exist('fbg', 'var')
fbg = 2.0                           % feedback gain to internal integrator
end

%
%   if there is an input soundfile specified, use it.  else, create a sin wave
%

if exist('inputFile', 'var')

fileSize = length(inputBuffer);

numSamples = 2.^(ceil(log2(fileSize(1))));  % round up to nearest power of 2

x = zeros(numSamples, 1);                   % zero pad if necessary

x(1:fileSize) = inputBuffer(:,1);           % if multi-channel, use left channel only

clear inputBuffer;                          % free this memory
clear fileSize;

t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time

else

if ~exist('numSamples', 'var')
numSamples = 65536                              % number of samples in simulation
end

if ~exist('Fs', 'var')
Fs = 44100                                      % (oversampled) sample rate
end

if ~exist('f0', 'var')
f0 = 261.6255653                                % input freq (middle C)
end

if ~exist('Amplitude', 'var')
Amplitude = 0.25                                % input amplitude
end

t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time
x = Amplitude*cos(2*pi*f0*t);                       % the input

end

sound(x, Fs);                                   % listen to input sound
pause;

y = zeros(1, numSamples);                       % the output (created and initialized for speed later)

if linearized_model
% artificial quantization noise for linearized model
% mean(q) = 0, var(q) = mean(q^2) = mean(y^2) - G^2*mean(v^2)
% does not have to be uniform or triangle p.d.f.
q = sqrt(6.0*(A^2 - G^2*mean_vv))*( rand(1, numSamples) - rand(1, numSamples) );
else
q = zeros(1, numSamples);
end

sum_yv = 0.0;
sum_vv = 0.0;

w = 0;
v = 0;
for n = 1:numSamples

if linearized_model

y(n) = G*v + q(n);                      % here the comparator is modelled as a little gain with additive noise

else

if (v >= 0)                             % the comparator
y(n) = +A;
else
y(n) = -A;
end

q(n) = y(n) - (sum_vv+1e-20)/(sum_yv+1e-20)*v;

end

sum_yv = sum_yv + y(n)*v;                   % collect some statistics on v
sum_vv = sum_vv +    v*v;

v = v + w  - fbg*y(n);                      % second integrator
w = w + x(n) - y(n);                        % first integrator

end

if ~linearized_model                            % don't recalculate this if using the linearized model
mean_yv = sum_yv/numSamples;
mean_vv = sum_vv/numSamples;
G = mean_yv/mean_vv;                        % the apparent comparator gain (assuming stationary input)
end

%
%
%
%     Y = ((G*z)*X + (z^3 - 2*z^2 + z)*Q) / (z^3 - 2*z^2 + (G*a+1)*z + G*(1-a))
%
%
%
Hx = freqz([0  0 G 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);
Hq = freqz([1 -2 1 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);

plot(t, y, 'b');
sound(y, Fs);                                   % this could sound pretty bad
pause;

Y = fft(fftshift(y .* kaiser(numSamples, 5.0)'));
Q = fft(fftshift(q .* kaiser(numSamples, 5.0)'));

f = linspace(0.0, (numSamples/2-1)/numSamples*Fs, numSamples/2);

plot(f, 20*log10(abs(Y(1:numSamples/2)) + 1e-10), 'b');
hold on;
plot(f, 20*log10(abs(Q(1:numSamples/2)) + 1e-10), 'r');
plot(f, 20*log10(abs(Hq) + 1e-10), 'g');
axis([0 Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Q(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hx(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 110]);
hold off;