Thursday, November 8, 2018

fft - Artifacts signal conversion from float to integer


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 figure with signal results



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



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()

Dithering flattens the noise floor
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()

Residual and its spectrum
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.


No comments:

Post a Comment

periodic trends - Comparing radii in lithium, beryllium, magnesium, aluminium and sodium ions

Apparently the of last four, $\ce{Mg^2+}$ is closest in radius to $\ce{Li+}$. Is this true, and if so, why would a whole larger shell ($\ce{...