English: Spectrum of a 750 Hz pure tone sampled at 48 kHz and quantized to 4 bits using triangular dither and noise shaping.
% quantization_shaped_dither.m
% Quantization demo using shaped dither.
% Authors: Fatbag
% License: Public domain (no warranties)
rand('state', 0xdeadbeef);
f_s = 48000; % sampling rate (samples/sec)
f0 = 750; % frequency of sine wave (Hz)
t_s = 1/f_s; % sampling period (sec/sample)
w0 = 2*pi*f0; % angular frequency of sine wave (Hz)
t = (0 : 4095) * t_s;
# x = round(16384*sin(w0*t) + rand(1, length(t)) - 0.5 + rand(1, length(t)) - 0.5);
x = zeros(1, length(t));
% FIR filter coefficients taken from Audacity (as of 2015-11-08)
h = [1, 2.033, -2.165, 1.959, -1.590, 0.6149];
prev_error = zeros(1, 5);
for i = 1 : length(t)
ideal = 15*sin(w0*t(i)) + dot(h(2:end), prev_error);
x(i) = round(ideal + rand() - 0.5 + rand() - 0.5);
prev_error = [ideal - x(i), prev_error(1:4)];
endfor
figure 1;
plot(t(1:ceil(2*f_s/f0)), x(1:ceil(2*f_s/f0)));
axis([0, ceil(2*f_s/f0)*t_s, -20, 20]);
grid on;
xlabel('t');
ylabel('x(t)');
figure 2;
X = fft(x(1:4096));
X = X(1:2048);
m = 20*log10(arrayfun(@norm, X));
m -= max(m);
m(~isfinite(m)) = -1000;
f = (0:2047)/2048*f_s/2;
plot(f, m);
axis([0, f_s/2, -100, 0]);
grid on;
xlabel('Frequency (Hz)');
ylabel('Relative power (dB)');
print 'quantization_shaped_dither_spectrum_uncropped.png';
figure 3;
stem(0 : (length(h)-1), h);
figure 4;
[H, w] = freqz(h);
m = arrayfun(@norm, H);
m = 20*log10(m);
plot(w, m);
axis([0, pi, -10, 50]);