t = linspace(0, (N-1)/fs, N);

y = Afs * cos(2*pi*fx*t);

figure('Units', 'normalized', 'Position', [0.05 0.15 0.85 0.7])

plot(t, y), title('Signal y'); axis tight

s = 1/Nfft * abs(fft(y, Nfft));

f = linspace(0, fs/2, Nfft/2+1);

sy = [s(1), 2*s(2:end-1), s(end)];

plot(f, pow2db(sy)), grid minor

title('One-sided power spectrum of y'), xlabel('Frequency - Hz');

yy = Afs * cos(2*pi*fx*t);

plot(t, yy), title('Signal yy'); axis tight

s = 1/Nfft * abs(fft(yy, Nfft));

syy = [s(1), 2*s(2:end-1), s(end)];

plot(f, pow2db(syy)), grid minor

title('One-sided power spectrum of yy'), xlabel('Frequency - Hz');

y1 = round(y/delta)*delta;

plot(t(1:n), y1(1:n)), title('Signal y1'); axis tight

ylabel('Amplitude'); grid minor

plot(t(1:n), y(1:n), '--')

legend({'Quantized', 'Original', 'Error'})

s = 1/Nfft * abs(fft(y1, Nfft));

sy1 = [s(1), 2*s(2:end-1), s(end)];

plot(f, pow2db(sy1)), grid minor

title('One-sided power spectrum of y1'), xlabel('Frequency - Hz');

pRMSyy = mag2db(rms(yy));

pRMSy1 = mag2db(rms(y1));

fprintf('Power of signal y = %g dBW\n', pRMSy)

Power of signal y = -3.01267 dBW

fprintf('Power of signal yy = %g dBW\n', pRMSyy)

Power of signal yy = -3.0103 dBW

fprintf('Power of signal y1 = %g dBW\n', pRMSy1)

Power of signal y1 = -2.22132 dBW