How would sigma-delta modulation work in software?

Question

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?

Answer 1

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.

Answer 2

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.)

Answer 3

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')

    [inputBuffer, Fs] = audioread(inputFile);

    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;

Answer 4

Apple mac app store HI-Player use delta sigma modulator and fir / iir filter upResample convert pcm to dsd signal