Time series dBFS plot output modification - current output plot not as expected (matplotlib)

Question

I'm trying to plot the Amplitude (dBFS) vs. Time (s) plot of an audio (.wav) file using matplotlib. I managed to do that with the following code:

def convert_to_decibel(sample):
    ref = 32768                    # Using a signed 16-bit PCM format wav file. So, 2^16 is the max. value.
    if sample!=0:
        return 20 * np.log10(abs(sample) / ref)

    else:
        return 20 * np.log10(0.000001)


from scipy.io.wavfile import read as readWav
from scipy.fftpack import fft

import matplotlib.pyplot as gplot1
import matplotlib.pyplot as gplot2
import numpy as np
import struct
import gc

wavfile1 = '/home/user01/audio/speech.wav'

wavsamplerate1, wavdata1 = readWav(wavfile1)
wavdlen1 = wavdata1.size
wavdtype1 = wavdata1.dtype

gplot1.rcParams['figure.figsize'] = [15, 5]
pltaxis1 = gplot1.gca()
gplot1.axhline(y=0, c="black")
gplot1.xticks(np.arange(0, 10, 0.5))
gplot1.yticks(np.arange(-200, 200, 5))
gplot1.grid(linestyle = '--')
wavdata3 = np.array([convert_to_decibel(i) for i in wavdata1], dtype=np.int16)
yvals3 = wavdata3
t3 = wavdata3.size / wavsamplerate1
xvals3 = np.linspace(0, t3, wavdata3.size)
pltaxis1.set_xlim([0, t3 + 2])
pltaxis1.set_title('Amplitude (dBFS) vs Time(s)')
pltaxis1.plot(xvals3, yvals3, '-')

which gives the following output:

I had also plotted the Power Spectral Density (PSD, in dBm) using the code below:

from scipy.signal import welch as psd            # Computes PSD using Welch's method.

fpsd, wPSD = psd(wavdata1, wavsamplerate1, nperseg=1024)

gplot2.rcParams['figure.figsize'] = [15, 5]

pltpsdm = gplot2.gca()
gplot2.axhline(y=0, c="black")
pltpsdm.plot(fpsd, 20*np.log10(wPSD))
gplot2.xticks(np.arange(0, 4000, 400))
gplot2.yticks(np.arange(-150, 160, 10))
pltpsdm.set_xlim([0, 4000])
pltpsdm.set_ylim([-150, 150])
gplot2.grid(linestyle = '--')

which gives the output as:

The second output above, using the Welch's method plots a more presentable output. The dBFS plot though informative is not very presentable IMO. Is this because of:

1. the difference in the domains (time in case of 1st output vs frequency in the 2nd output)?
2. the way plot function is implemented in pyplot?

Also, is there a way I can plot my dBFS output as a peak-to-peak style of plot just like in my PSD (dBm) plot rather than a dense stem plot?

Would be much helpful and would appreciate any pointers, answers or suggestions from experts here as I'm just a beginner with matplotlib and plots in python in general.

Answer 1

TLNR

• This has nothing to do with pyplot.
• The frequency domain is different from the time domain, but that's not why you didn't get what you want.
• The calculation of dbFS in your code is wrong.

You should frame your data, calculate RMSs or peaks in every frame, and then convert that value to dbFS instead of applying this transformation to every sample point.

---

When we talk about the amplitude, we are talking about a periodic signal. And when we read in a series of data from a sound file, we read in a series of sample points of a signal(may be or be not periodic). The value of every sample point represents a, say, voltage value, or sound pressure value sampled at a specific time.

We assume that, within a very short time interval, maybe 10ms for example, the signal is stationary. Every such interval is called a frame.

Some specific function is applied to each frame usually, to reduce the sudden change at the edge of this frame, and these functions are called window functions. If you did nothing to every frame, you added rectangle windows to them.

An example: when the sampling frequency of your sound is 44100Hz, in a 10ms-long frame, there are 44100*0.01=441 sample points. That's what the nperseg argument means in your psd function but it has nothing to do with dbFS.

Given the knowledge above, now we can talk about the amplitude.

There are two methods a get the value of amplitude in every frame:

• The most straightforward one is to get the maximum(peak) values in every frame.
• Another one is to calculate the RMS(Root Mean Sqaure) of every frame.

After that, the peak values or RMS values can be converted to dbFS values.

Let's start coding:

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

# Determine full scall(maximum possible amplitude) by bit depth
bit_depth = 16
full_scale = 2 ** bit_depth

# dbFS function
to_dbFS = lambda x: 20 * np.log10(x / full_scale)

# Read in the wave file
fname = "01.wav"
fs,data = wavfile.read(fname)

# Determine frame length(number of sample points in a frame) and total frame numbers by window length(how long is a frame in seconds)
window_length = 0.01 
signal_length = data.shape[0]
frame_length = int(window_length * fs)
nframes = signal_length // frame_length

# Get frames by broadcast. No overlaps are used.
idx = frame_length * np.arange(nframes)[:,None] + np.arange(frame_length)
frames = data[idx].astype("int64") # Convert to in 64 to avoid integer overflow

# Get RMS and peaks
rms = ((frames**2).sum(axis=1)/frame_length)**.5
peaks = np.abs(frames).max(axis=1)

# Convert them to dbfs
dbfs_rms = to_dbFS(rms)
dbfs_peak = to_dbFS(peaks)

# Let's start to plot

# Get time arrays of every sample point and ever frame
frame_time = np.arange(nframes) * window_length
data_time = np.linspace(0,signal_length/fs,signal_length)

# Plot
f,ax = plt.subplots()
ax.plot(data_time,data,color="k",alpha=.3)

# Plot the dbfs values on a twin x Axes since the y limits are not comparable between data values and dbfs
tax = ax.twinx()
tax.plot(frame_time,dbfs_rms,label="RMS")
tax.plot(frame_time,dbfs_peak,label="Peak")
tax.legend()
f.tight_layout()

# Save serval details
f.savefig("whole.png",dpi=300)
ax.set_xlim(1,2)
f.savefig("1-2sec.png",dpi=300)
ax.set_xlim(1.295,1.325)
f.savefig("1.2-1.3sec.png",dpi=300)

The whole time span looks like(the unit of the right axis is dbFS):

And the voiced part looks like:

You can see that the dbFS values become greater while the amplitudes become greater at the vowel start point:

Answer 2

I don't necassarily think there is anything wrong with your first plot, but what you might want to do is reduce the noise. Fourier transform your signal to the frequency domain, remove the higher frequencies with some threshold and inverse fourier transform back to the time domain. Try a couple different thresholds and see what you like. Here's a blogg post that should contain sufficent code examples.

https://towardsdatascience.com/noise-cancellation-with-python-and-fourier-transform-97303314aa71