时间序列 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 上并不是很像样。这是因为:
- 域的差异(第一个输出的时间与第二个输出的频率)?
- 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 值变大,而振幅在元音起始点变大:
我正在尝试使用 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 上并不是很像样。这是因为:
- 域的差异(第一个输出的时间与第二个输出的频率)?
- 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 值变大,而振幅在元音起始点变大: