close all; clear; clc

%% Parameter

w = 2;
fs = 44100;
N = 3001;

Q = 2/2^w;

%% Eingangssignale

f = 1000;
A = 0.5;

x = A * sin(2*pi*[0:N-1]*f/fs);


%% Dither erzeugen

% d(n) = d(n) - d(n-1)

d_rect = Q/2 * (2*rand(1,N+1)-1);
d_HP = d_rect(2:N+1) - d_rect(1:N);

% hist(d_rect)
% hist(d_HP)
% max(d_HP)

%% Quantisierung auf w bit

bUseDither = true;

for iC = 1:N

    if bUseDither
        xqr(iC) = Q * floor( (x(iC)+d_HP(iC)) / Q + 0.5 );
    else
        xqr(iC) = Q * floor(x(iC) / Q + 0.5 );
    end
end

er = x -xqr;

%% Zeigsignal darstellen

figure('Name', 'Time signal');
subplot(211)
plot(0:N-1, x, 'b.-', 'DisplayName', 'in')
grid on; hold on;
plot(0:N-1, xqr, 'ro-', 'DisplayName', 'quantized')
legend('show', 'Location', 'Best');
xlabel('n \rightarrow')

subplot(212)
plot(0:N-1, er, 'k', 'DisplayName', 'Quantisation Error')

%% Spektrum

Nfft = 1024;

Xdb = 10*log10( 1/(Nfft) * abs( fft(x,Nfft) ) );
Xqrdb = 10*log10( 1/(Nfft) * abs( fft(xqr,Nfft) ) );

fFFT = fs * (0:Nfft-1)/Nfft;

figure('Name', 'Spectrum');
plot(fFFT,Xdb, 'b.-', 'DisplayName', 'in')
grid on; hold on;
plot(fFFT,Xqrdb, 'ro-', 'DisplayName', 'in')
legend('show', 'Location', 'Best');
xlabel('n \rightarrow')