1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Widrow, B. came out with a book on amplitude about 10 years ago. He also published some papers on the subject. you can find the references off his web page www-isl.stanford.edu/~widrow $\endgroup$ – user28715 Feb 5 '18 at 19:01
  • $\begingroup$ I suggest bugfixes to your script: Change x = np.arange(l) / l to x = np.arange(l) / float(l), change each 10 * np.log10 to 20 * np.log10, add to the end fig.show(). $\endgroup$ – Olli Niemitalo Feb 14 '18 at 9:03
  • $\begingroup$ Just as a note, @OlliNiemitalo and I have spotted a couple more things you might want to revise in your code / thinking in our responses below. $\endgroup$ – A_A Feb 14 '18 at 12:32
6
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Thank you. I did the same. I do not mind "being wrong" and I would not like this to be perceived as an argument. I would however like to understand what is going on here. $\endgroup$ – A_A Feb 14 '18 at 11:24
  • $\begingroup$ Thank you. +1, indeed, that is a key detail. I have updated my answer as well. $\endgroup$ – A_A Feb 14 '18 at 12:05
3
$\begingroup$

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.

Time Frequency 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.

Time Freq plot of the quantised 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)));

Time freq plot slightly quantised

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

Time freq plot, heavy distortion

...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:

Residuals

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:

Residuals 2

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:

Residuals 3

$\endgroup$
  • $\begingroup$ It's not clipping. (sine / sine.max() * max_int).min() and (sine / sine.max() * max_int).max() return -32767.0 and 32767.0, which both fit into a signed 16-bit integer. $\endgroup$ – Olli Niemitalo Feb 14 '18 at 9:13
  • $\begingroup$ @OlliNiemitalo thanks for letting me know. In my attempt I tried to replicate the error. Quantisation noise doesn't have such "well behaved" spectrum. The harmonics that we see in the question suggest a specific waveform and are very regular. The spectrum of the quantisation residual is not that regular. In fact, plain conversion did not produce any artifacts, it was only after clipping was introduced that the artefacts appeared... $\endgroup$ – A_A Feb 14 '18 at 9:45
  • $\begingroup$ @OlliNiemitalo ...Therefore, what I would like to suggest is to try and replicate the error, with and without dithering so that it is clear what is causing what. I am not saying that dithering the waveform is not advisable here of course. Rather, I think that it would be a very good additional step to take anyway. $\endgroup$ – A_A Feb 14 '18 at 9:51
  • $\begingroup$ I added stuff to my answer to show how the well behaved spectrum comes from the quantization error. $\endgroup$ – Olli Niemitalo Feb 14 '18 at 10:39
  • 1
    $\begingroup$ I spotted the difference: Jan-Bert quantizes before windowing. $\endgroup$ – Olli Niemitalo Feb 14 '18 at 11:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.