2
$\begingroup$

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?

$\endgroup$
2
$\begingroup$
  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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Is the dithering step before modulating the signal with PDM? $\endgroup$ – Fabio K Nov 5 '13 at 15:05
  • $\begingroup$ 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. $\endgroup$ – adiblol Nov 5 '13 at 15:22
0
$\begingroup$

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

Block Diagram

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

Spectrum

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – robert bristow-johnson May 2 at 1:12
  • $\begingroup$ 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 $\endgroup$ – Dan Boschen May 2 at 1:15
0
$\begingroup$

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('A', '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;
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.