% 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]);