python 中的频谱图使用 numpy
Spectrogram in python using numpy
我需要使用 numpy 制作频谱图。我将 1 秒的音频分成 0.02 秒的块。然后我使用 numpy 计算 FFT 并将其重新组合成一个图像。结果很差。
这是使用 matplotlib specgram 函数生成的频谱图：
这是我的 'spectrogram':
这是我的代码：
spect_frags = []
transform = []
for x in range(0, 8000, 160):
spect_frags.append(spect_sample[x:x + 160])
for sample in spect_frags:
transform.append(abs(np.fft.fft(sample).real)[0:np.fft.fft(sample).real.size//4])
我剪掉了 3/4 的频率，因为我现在不需要它们。
我不知道为什么分辨率会有如此大的差异。我该如何改进它？
MCVE 频谱图
您可以使用以下代码重新创建 specgram 的粗略估计：
import numpy as np
from scipy.io import wavfile
from scipy import fft
import matplotlib.pyplot as plt
# Read some sample file (replace with your data):
rate, data = wavfile.read('./data/aaaah.wav')
# rate=48000, data.shape=(46447, 2) ~ almost 1s of stereo signal
# Spectrogram estimation:
N = 256
S = []
for k in range(0, data.shape[0]+1, N):
x = fft.fftshift(fft.fft(data[k:k+N,0], n=N))[N//2:N]
# assert np.allclose(np.imag(x*np.conj(x)), 0)
Pxx = 10*np.log10(np.real(x*np.conj(x)))
S.append(Pxx)
S = np.array(S)
# Frequencies:
f = fft.fftshift(fft.fftfreq(N, d=1/rate))[N//2:N]
# array([ 0. , 187.5, 375. , ..., 23625. , 23812.5])
# Spectrogram rendering:
plt.imshow(S.T, origin='lower')
它输出：
当 specgram 渲染时：
_ = plt.specgram(data[:,0])
此 MCVE 与 specgram 不同，因为应缩放轴以正确反映时间和频率，并且没有移动 windowing。更准确地说：
- x-axis表示长度为
N=256的时间块索引；
- y-axis是正半平面FFT指数(
N//2=128)，注意fftshift到assemblefft后的频谱的使用；
- 实际频率可以使用采样率和
fftfreq，在 specgram 中，它的范围从 0 到 1，因为这种方法不一定知道信号采样率；
- 没有 window 重叠（使用了独立的连续块），这就是为什么 MCVE 不如
specgram. 平滑的原因
功率估计
另请注意，取复数的实部与取大小不同。主要是，当你写：
abs(np.fft.fft(sample).real)
您没有采用复数的范数，但由于 .real 调用，您完全删除了复数部分。
你应该estimate the power using product of conjugates:
10*np.log10(np.real(x*np.conj(x)))
然后使用abs将complex类型（或者只保留real部分，因为复数部分必须为空）变成float。最后，您可以使用十进制对数缩放 Decibel。
完整性检查
可以查看FFT的结果确实是复数类型，有显着的复数部分（去掉会丢失信息）：
x
# array([-1.56000000e+02-0.00000000e+00j, -3.94271344e+01+1.17935735e+02j,
# 4.03754070e+01+4.14695163e+01j, 1.71510716e+01+1.26920718e+01j,
# 2.15523795e+01-2.07362424e+00j, -3.03847433e+00-1.22767815e+01j,
# -4.56347533e+00-7.36380957e-01j, -1.28048283e+01-6.80931256e+00j,
# -2.22781473e+01+1.12096897e+01j, -1.13788549e+01+2.54314337e+01j,
# ...])
并且共轭的乘积确实有一个空复数部分（但仍然是 complex 类型）：
x*np.conj(x)
# array([2.43360000e+04+0.j, 1.54633365e+04+0.j, 3.34989427e+03+0.j,
# 4.55247945e+02+0.j, 4.68804979e+02+0.j, 1.59951690e+02+0.j,
# 2.13675640e+01+0.j, 2.10330365e+02+0.j, 6.21972990e+02+0.j,
# 7.76236159e+02+0.j, 1.05846430e+03+0.j, 6.54663598e+02+0.j,
# 6.95792718e+01+0.j, 6.03013130e+01+0.j, 1.11620428e+01+0.j,
# ...])
您可以通过断言以下内容来确保这始终为真（完整性检查）：
assert np.allclose(np.imag(x*np.conj(x)), 0)
为了有效地使用 JAX I found it useful to adapt @jlandercy 的解决方案来避免显式循环并添加一些简单的 Hann 窗口。我在此过程中放弃了对立体声输入的支持，尽管我相当确定计算 STFT 的全部意义在于傅里叶变换的非线性意味着我在执行此操作之前通过将立体声信号折叠为单声道来破坏一些信息。
我也没有扩展他们的答案以包括重叠 windows 以保持其简单性和简洁性。可能有一些 numpy 女巫知道如何高效地 np.concatenate 拼接在一起或从 wins 或类似的最内层维度插入非重叠信号样本的切片。遗憾的是，我本人并没有这方面的黑魔法知识，所以功能被省略了。
import matplotlib.pyplot as plt
import jax.numpy as jnp
def stft(a, n_fft=128, window=jnp.hanning):
n = n_fft
rpad = n - a.shape[-1] % n
wins = jnp.pad(a, (0, rpad)).reshape(-1, n) * window(n)
fftc = jnp.fft.fftshift(jnp.fft.fft(wins, n=n))[..., n // 2 : n]
fftr = jnp.real(fftc * jnp.conj(fftc))
return fftr
audio = # buffer single-channel floating-point samples from somewhere...
_ = plt.imshow(stft(audio, 1024).T[:, -512:], cmap="viridis")
我需要使用 numpy 制作频谱图。我将 1 秒的音频分成 0.02 秒的块。然后我使用 numpy 计算 FFT 并将其重新组合成一个图像。结果很差。
这是使用 matplotlib specgram 函数生成的频谱图：
这是我的 'spectrogram':
这是我的代码：
spect_frags = []
transform = []
for x in range(0, 8000, 160):
spect_frags.append(spect_sample[x:x + 160])
for sample in spect_frags:
transform.append(abs(np.fft.fft(sample).real)[0:np.fft.fft(sample).real.size//4])
我剪掉了 3/4 的频率，因为我现在不需要它们。 我不知道为什么分辨率会有如此大的差异。我该如何改进它？
MCVE 频谱图
您可以使用以下代码重新创建 specgram 的粗略估计：
import numpy as np
from scipy.io import wavfile
from scipy import fft
import matplotlib.pyplot as plt
# Read some sample file (replace with your data):
rate, data = wavfile.read('./data/aaaah.wav')
# rate=48000, data.shape=(46447, 2) ~ almost 1s of stereo signal
# Spectrogram estimation:
N = 256
S = []
for k in range(0, data.shape[0]+1, N):
x = fft.fftshift(fft.fft(data[k:k+N,0], n=N))[N//2:N]
# assert np.allclose(np.imag(x*np.conj(x)), 0)
Pxx = 10*np.log10(np.real(x*np.conj(x)))
S.append(Pxx)
S = np.array(S)
# Frequencies:
f = fft.fftshift(fft.fftfreq(N, d=1/rate))[N//2:N]
# array([ 0. , 187.5, 375. , ..., 23625. , 23812.5])
# Spectrogram rendering:
plt.imshow(S.T, origin='lower')
它输出：
当 specgram 渲染时：
_ = plt.specgram(data[:,0])
此 MCVE 与 specgram 不同，因为应缩放轴以正确反映时间和频率，并且没有移动 windowing。更准确地说：
- x-axis表示长度为
N=256的时间块索引； - y-axis是正半平面FFT指数(
N//2=128)，注意fftshift到assemblefft后的频谱的使用； - 实际频率可以使用采样率和
fftfreq，在specgram中，它的范围从 0 到 1，因为这种方法不一定知道信号采样率； - 没有 window 重叠（使用了独立的连续块），这就是为什么 MCVE 不如
specgram. 平滑的原因
功率估计
另请注意，取复数的实部与取大小不同。主要是，当你写：
abs(np.fft.fft(sample).real)
您没有采用复数的范数，但由于 .real 调用，您完全删除了复数部分。
你应该estimate the power using product of conjugates:
10*np.log10(np.real(x*np.conj(x)))
然后使用abs将complex类型（或者只保留real部分，因为复数部分必须为空）变成float。最后，您可以使用十进制对数缩放 Decibel。
完整性检查
可以查看FFT的结果确实是复数类型，有显着的复数部分（去掉会丢失信息）：
x
# array([-1.56000000e+02-0.00000000e+00j, -3.94271344e+01+1.17935735e+02j,
# 4.03754070e+01+4.14695163e+01j, 1.71510716e+01+1.26920718e+01j,
# 2.15523795e+01-2.07362424e+00j, -3.03847433e+00-1.22767815e+01j,
# -4.56347533e+00-7.36380957e-01j, -1.28048283e+01-6.80931256e+00j,
# -2.22781473e+01+1.12096897e+01j, -1.13788549e+01+2.54314337e+01j,
# ...])
并且共轭的乘积确实有一个空复数部分（但仍然是 complex 类型）：
x*np.conj(x)
# array([2.43360000e+04+0.j, 1.54633365e+04+0.j, 3.34989427e+03+0.j,
# 4.55247945e+02+0.j, 4.68804979e+02+0.j, 1.59951690e+02+0.j,
# 2.13675640e+01+0.j, 2.10330365e+02+0.j, 6.21972990e+02+0.j,
# 7.76236159e+02+0.j, 1.05846430e+03+0.j, 6.54663598e+02+0.j,
# 6.95792718e+01+0.j, 6.03013130e+01+0.j, 1.11620428e+01+0.j,
# ...])
您可以通过断言以下内容来确保这始终为真（完整性检查）：
assert np.allclose(np.imag(x*np.conj(x)), 0)
为了有效地使用 JAX I found it useful to adapt @jlandercy 的解决方案来避免显式循环并添加一些简单的 Hann 窗口。我在此过程中放弃了对立体声输入的支持，尽管我相当确定计算 STFT 的全部意义在于傅里叶变换的非线性意味着我在执行此操作之前通过将立体声信号折叠为单声道来破坏一些信息。
我也没有扩展他们的答案以包括重叠 windows 以保持其简单性和简洁性。可能有一些 numpy 女巫知道如何高效地 np.concatenate 拼接在一起或从 wins 或类似的最内层维度插入非重叠信号样本的切片。遗憾的是，我本人并没有这方面的黑魔法知识，所以功能被省略了。
import matplotlib.pyplot as plt
import jax.numpy as jnp
def stft(a, n_fft=128, window=jnp.hanning):
n = n_fft
rpad = n - a.shape[-1] % n
wins = jnp.pad(a, (0, rpad)).reshape(-1, n) * window(n)
fftc = jnp.fft.fftshift(jnp.fft.fft(wins, n=n))[..., n // 2 : n]
fftr = jnp.real(fftc * jnp.conj(fftc))
return fftr
audio = # buffer single-channel floating-point samples from somewhere...
_ = plt.imshow(stft(audio, 1024).T[:, -512:], cmap="viridis")