时间序列 dBFS 图输出修改 - 当前输出图不符合预期(matplotlib)

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

我正在尝试使用 matplotlib 绘制音频 (.wav) 文件的 Amplitude (dBFS) vs. Time (s) 图。我设法使用以下代码做到了这一点:

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, '-')

给出以下输出:

我还使用以下代码绘制了 Power Spectral Density (PSD, in dBm)

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 = '--')

输出为:

上面的第二个输出,使用 Welch 的方法绘制了一个更像样的输出。 dBFS 图虽然提供了丰富的信息,但在 IMO 上并不是很像样。这是因为:

  1. 域的差异(第一个输出的时间与第二个输出的频率)?
  2. pyplot中plot函数的实现方式?

此外,有没有一种方法可以将我的 dBFS 输出绘制为 峰对峰样式的图 就像我的 PSD (dBm) 情节而不是密集的 干图 ?

会非常有帮助,并且会感谢这里专家的任何指示、答案或建议,因为我只是 matplotlib 的初学者,一般来说 python 中的绘图。

我不一定认为你的第一个情节有什么问题,但你可能想要做的是减少噪音。将信号傅里叶变换到频域,用一些阈值去除较高的频率,然后傅里叶逆变换回时域。尝试几个不同的阈值,看看你喜欢什么。这是一个 blogg post,应该包含足够的代码示例。

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

TLNR

  • 这与pyplot无关。
  • 频域与时域不同,但这不是您没有得到想要的结果的原因。
  • 你代码中的dbFS计算错误。

您应该构建您的数据,计算每一帧中的 RMS 或峰值,然后将该值转换为 dbFS,而不是将此转换应用于每个样本点。


当我们谈论幅度时,我们谈论的是周期信号。当我们从声音文件中读取一系列数据时,我们读取了信号的一系列 样本点 (可能是周期性的,也可能不是周期性的)。每个采样点的值代表一个电压值,或者在特定时间采样的声压值。

我们假设,在很短的时间间隔内,例如 10ms,信号是静止的。每一个这样的间隔称为一个.

通常对每一帧应用一些特定的函数,以减少这一帧边缘的突然变化,这些函数称为window函数。如果你对每一帧什么都不做,你就给它们添加了矩形windows。

举个例子:当你的声音采样频率为44100Hz时,在一个10ms长的帧中,有44100*0.01=441个采样点。这就是 nperseg 参数在 psd 函数中的含义,但它与 dbFS 无关。

有了上面的知识,现在我们可以谈谈振幅了。

获取每一帧振幅值的方法有两种:

  • 最直接的方法是获取每一帧中的最大(峰值)值。
  • 另一种是计算每一帧的RMS(Root Mean Sqaure)。

之后,可以将峰值或RMS值转换为dbFS值。

让我们开始编码:

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)

整个时间跨度是这样的(右轴的单位是dbFS):

有声部分看起来像:

您可以看到 dbFS 值变大,而振幅在元音起始点变大: