Artifacts signal conversion from float to integer

Question

I have a script where i want to change the noise floor to about 96 dB (visual purpose). The script containing basically a sine wave with after wart do FFT over it.

The script is in python and as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / l
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
sine_96 = (sine / sine.max() * max_int).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(10 * np.log10(S[:l // 2]))
ax[1, 1].plot(10 * np.log10(S_96[:l // 2]))

The weird thing is when i change the signal to 16 bit integer (or 16 bit integer and back to floats again) i get some artefacts what i don't understand. See also

left are the floating point functions where i get about an 150 dB noise floor. and on the right the version where the signal is changed to 16 bit integer. The result should be the same signal but with a noise floor around 96 dB.

As i can see it the noise floor is getting higher but so do the peak. I also get's repetitive peaks. I have two question how can i avoid these peaks and what is the background of this?

Answer 1

The quantization error is dependent on the input signal. In your case the input signal is a sine wave and the quantization error is a combination of its 8 odd harmonics that fit into the available bandwidth. You can reduce the dependency of the quantization error to the input signal by adding triangular dither before quantizing: Add to each sample a uniformly distributed random number between 0 and the quantization step, and subtract from the result another uniformly distributed random number between 0 and the quantization step. Then quantize. This way you'll get a flat spectral noise floor (no distortion of the signal) and no modulation of the quantization noise by the signal.

For further reading on dither: Robert Alexander Wannamaker's PhD thesis "The Theory of Dithered Quantization", 2003.

Modifying your script to include dithering fixes the problem:

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / float(l)
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
dither = [random.triangular(-1, 1, 0) for _ in xrange(l)]
sine_96 = (sine / sine.max() * (max_int - 1) + dither).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(20 * np.log10(S[:l // 2]))
ax[1, 1].plot(20 * np.log10(S_96[:l // 2]))
fig.show()

Figure 1. Top left: windowed input signal, top right: dithered, quantized, and windowed input signal, bottom left: frequency spectrum of the windowed input signal (vertical axis in dB scale), bottom right: frequency spectrum of the dithered, quantized, and windowed input signal (vertical axis in dB scale).

Let's plot the quantization residual (error) of the undithered case before windowing, and the frequency spectrum of the residual after windowing it:

sine_96 = (sine / sine.max() * max_int).astype(np.int16)
residual = sine_96 - (sine / sine.max() * max_int)
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
ax[0].plot(residual)
R = np.abs(fft(residual * w))
ax[1].plot(20 * np.log10(R[:l // 2]))
fig.show()

Figure 2. Quantization residual and its spectrum (window function applied before transforming).

Quantization is equivalent to a memoryless waveshaper shaped like a staircase. For a periodic input signal a memoryless waveshaper can only produce harmonics of the input signal, because the resulting output signal will be identically periodic:

$$f(x + s) = f(x) \text{ for all } x \in \mathbf{R}\\ \Rightarrow g\big(f(x + s)\big) = g\big(f(x)\big) \text{ for all } x \in \mathbf{R}$$

Here the input signal is a sine wave with a period of 128/4 = 32 samples. If the period would not be an integer number of samples, then also other frequencies than the harmonics would be generated, because quantization operates in discrete time and non-integer periodicities don't exist in discrete time.

Answer 2

Your quantised signal is overflowing the range of a 16bit integer and it is clipping.

I will replicate and explain the situation here, in Octave. Any inferences are easily transferable to Python.

Fs = 2000; %Sampling frequency in Hz
T = 2; % Total duration of the a signal in seconds
t = 0:(1./Fs):(T - (1./Fs)); %Time vector in steps of Ts = 1./Fs;
p = 2.0*pi.*t; %Phase vector in radians.
%
% Define the signal
%
s = hamming(length(t))'.*sin(50.*p); %Hamming modulated, 50Hz signal
%
% Let's have a quick look at it
%
subplot(121);plot(s);xlabel('Time');ylabel('Amplitude');grid on;subplot(122);semilogy(abs(fft(s)));xlabel('Frequency bin');ylabel('Amplitude');grid on;

This produces something predictable and inline with what you present on the left of your plot.

The only difference to this is that I am using semilogy for the Discrete Fourier Transform which is a better option to plotting spectra than a simple linear plot because of the wild difference in sums that can possibly be produced by the DFT.

So, no surprises here, 50Hz clean as a whistle. Now, let's quantise s in a straightforward way.

q = int16(round(s.*((2^15)-1))); % Signed integer max is 2^15 MINUS 1

Now, substitute q for s in our little plotting line above and let's have another look at the signal.

Still a whistle but not so clean now. And in fact, if you were to push it slightly overboard the int16 range, you get:

q = int16(round(s.*((2^15)+20)));

Although, to the human eye it is not yet perceptible, the waveform is clipping as it is evident from the emergence of those odd harmonics.

And if you keep pushing...

q = int16(round(s*((2^15)+8000)));

...results are now almost identical.

The clipping is the result of the intXX functions. If you tried to convert to rounded integers "manually", this would probably show up as wrap-around error.

Hope this helps.

EDIT:

The additional work performed by Olli Niemitallo further below is enough to convince me to look more into this. Still in Octave, I decided to get a plot of the quantisation residual. The point that I perceive as surprising is the magnitude of the quantisation noise at 16bit. Here is the result that I get in Octave:

Fs = 2000; %Sampling frequency in Hz
T = 2; % Total duration of the a signal in seconds
t = 0:(1./Fs):(T - (1./Fs)); %Time vector in steps of Ts = 1./Fs;
p = 2.0*pi.*t; %Phase vector in radians.
%
% Define the signal
%
s = hamming(length(t))'.*sin(50.*p); %Hamming modulated, 50Hz signal
%
%Up to here the script is the same and it produces the same output (of course). In the following two lines, I quantise in a very straightforward way and obtain the difference between the non-quantised and quantised waveforms
%
sq = int16(round(s*((2^15)-1))); %Signal Quantised
sqdq = double(sq)./((2^15)-1); % Signal quantised and then "de-quantised", notice the division by the absolute maximum, the maximum of the frame of reference, not the maximum of the waveform, this is important.
residual = sqdq-s;

Now, if I take a look at the residual in both domains, it looks like this:

Is there "structure" in this spectrum? Yes there is (familiar) structure in this spectrum, no doubt, but there is a "structure" at an order of magnitude of $e^{-5}$ and the original post shows considerable spikes.

EDIT 2:

Olli Niemitalo has spotted what is making the difference. Indeed, if we window the quantised waveform, then we see things differently:

Fs = 2000; %Sampling frequency in Hz
T = 2; % Total duration of the a signal in seconds
t = 0:(1./Fs):(T - (1./Fs)); %Time vector in steps of Ts = 1./Fs;
p = 2.0*pi.*t; %Phase vector in radians.
%
% Define the signal
%
s = sin(50.*p); %Hamming modulated, 50Hz signal
%
%Up to here the script is the same and it produces the same output (of course). In the following two lines, I quantise in a very straightforward way and obtain the difference between the non-quantised and quantised waveforms
%
sq = int16(round(s*((2^15)-1))); %Signal Quantised (ONLY THE SINUSOID)
sq2 = int16(round(sq.*hamming(length(sq))')); %Quantised signal windowed
sq2dsq = double(sq2)/((2^15)-1); %"De-quantise"
residual = sq2dsq-(s.*hamming(length(s))');

Now, this paints a slightly different picture:

In this case, I think both myself and Olli Niemitallo were thrown off by what we see.

The real reason for the "harmonics" seems to be "modulation".

When you quantise, you do indeed wave shaping with a "jagged line". But when you window (modulate) after that, then the steps of that "jagged line" are not horizontal anymore but slanted.

And that could really be the reason behind those harmonics that we see. Neither clipping, nor quantisation, both of which are real concerns of course.

The only reservation that I still have is regarding the magnitude that this happens at but I do not think that there is any doubt regarding what is causing what we see.

EDIT 3:

Just adding the image from the comments which clarifies further the exact kind of processing and its impact on quantisation: