python - python中使用numpy的频谱图
问题描述
我需要使用 numpy 制作频谱图。我取 1 秒的音频并将其分成 0.02 秒的块。然后我使用 numpy 计算 FFT 并将其重新组合成一个图像。结果很差。
这是使用 matplotlib specgram 函数生成的频谱图:
这是我的“频谱图”:
这是我的代码:
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
因为轴应该被缩放以正确反映时间和频率,并且没有移动窗口。更确切地说:
- x轴代表长度的时间块索引
N=256
; - y轴为正半平面
N//2=128
FFTfftshift
指数fft
( - 使用采样率可以获得真实频率,并且
fftfreq
在specgram
它的范围内从 0 到 1,因为这种方法不一定知道信号采样率; - 没有窗口重叠(使用独立的连续块),这就是为什么 MCVE 比
specgram
.
功率估算
还要注意取复数的实部与取幅值不同。主要是,当你写:
abs(np.fft.fft(sample).real)
您没有采用复数的规范,但是由于.real
调用,您完全删除了复杂的部分。
10*np.log10(np.real(x*np.conj(x)))
然后使用abs
将complex
类型(或仅保留real
部分,因为复杂部分必须为空)转换为float
. 最后,您可以使用十进制对数来缩放分贝。
完整性检查
您可以检查 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)
推荐阅读
- azure-sql-database - Azure SQL 数据库用户未继承角色权限
- c# - LINQ 查询表达式是否优化?
- python - 如何使用python将完成的抓取扩展到第一页以上
- r - R - 预测 - 根据上一年预测每月值
- python - 有没有办法阻止 python turicreate 使用子目录作为训练数据?
- javascript - 使用 javascript 设置活动类的简单方法
- java - 使用 AES 256 和 SHA-2 加密
- devexpress - 命名空间“DevExpress.Utils.About”中不存在类型或命名空间名称“frmAbout”
- php - 警告:mysqli::query(): 无法获取 mysqli - PHP7.1 & MySQL5.7
- matlab - MATLAB 在 t = 0、0.001、0.01、1 和 10 处绘制一个图