如何从 Python 中的 FFT 得到 time/freq
How to get time/freq from FFT in Python
我在管理 FFT 数据时遇到了一点问题。我一直在寻找许多如何进行 FFT 的示例，但我无法从其中任何一个中得到我想要的东西。我有一个 44kHz 采样率的随机波形文件，我想每 X 毫秒获得 N 次谐波的幅度，假设 100 毫秒就足够了。我试过这段代码：
import scipy.io.wavfile as wavfile
import numpy as np
import pylab as pl
rate, data = wavfile.read("sound.wav")
t = np.arange(len(data[:,0]))*1.0/rate
p = 20*np.log10(np.abs(np.fft.rfft(data[:2048, 0])))
f = np.linspace(0, rate/2.0, len(p))
pl.plot(f, p)
pl.xlabel("Frequency(Hz)")
pl.ylabel("Power(dB)")
pl.show()
这是我使用的最后一个示例，我在 Whosebug 的某个地方找到了它。问题是，这得到了我想要的幅度，得到了频率，但根本没有时间。据我所知，FFT 分析是 3D 的，这是所有谐波的 "merged" 结果。我明白了：
X-axis = Frequency, Y-axis = Magnitude, Z-axis = Time (invisible)
根据我对代码的理解，t 是时间 - 看起来是这样，但代码中不需要 - 不过我们可能会需要它。 p 是幂（或幅度）数组，但它似乎是每个频率 f 的所有幅度的某个平均值，它是频率数组。我不想要 average/merged 值，我想要每 X 毫秒 N 次谐波的幅度。
长话短说，我们可以得到：所有频率的1级。
我们想要：N 个频率的所有幅度，包括出现特定幅度时的时间。
结果应类似于此数组：[时间、频率、振幅]
所以最后如果我们想要 3 个谐波，它看起来像：
[0,100,2.85489] #100Hz harmonic has 2.85489 amplitude on 0ms
[0,200,1.15695] #200Hz ...
[0,300,3.12215]
[100,100,1.22248] #100Hz harmonic has 1.22248 amplitude on 100ms
[100,200,1.58758]
[100,300,2.57578]
[200,100,5.16574]
[200,200,3.15267]
[200,300,0.89987]
不需要可视化，结果应该只是上面列出的数组（或 hashes/dictionaries）。
您似乎正在尝试实现 spectrogram，这是一系列功率谱估计，通常通过一系列（通常重叠的）FFT 来实现。由于您只有一个 FFT（频谱），因此您还没有时间维度。将 FFT 代码放入循环中，每次迭代处理一个样本块（例如 1024），连续块之间有 50% 的重叠。生成的频谱序列将是时间 v 频率 v 幅度的 3D 数组。
我不是 Python 人，但我可以给你一些伪代码，应该足以让你编码：
N = length of data input
N_FFT = no of samples per block (== FFT size, e.g. 1024)
i = 0 ;; i = index of spectrum within 3D output array
for block_start = 0 to N - block_start
block_end = block_start + N_FFT
get samples from block_start .. block_end
apply window function to block (e.g. Hamming)
apply FFT to windowed block
calculate magnitude spectrum (20 * log10( re*re + im*im ))
store spectrum in output array at index i
block_start += N_FFT / 2 ;; NB: 50% overlap
i++
end
根据@Paul R 的回答，scipy.signal.spectrogram 是 spectrogram function in scipy's signal processing module。
上面link的例子如下：
from scipy import signal
import matplotlib.pyplot as plt
# Generate a test signal, a 2 Vrms sine wave whose frequency linearly
# changes with time from 1kHz to 2kHz, corrupted by 0.001 V**2/Hz of
# white noise sampled at 10 kHz.
fs = 10e3
N = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
freq = np.linspace(1e3, 2e3, N)
x = amp * np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
#Compute and plot the spectrogram.
f, t, Sxx = signal.spectrogram(x, fs)
plt.pcolormesh(t, f, Sxx)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()
编辑： 哦，看起来这个 returns 值，但它们根本不适合音频文件。尽管它们可以用作频谱图上的幅度，但它们在您可以在许多音乐播放器中看到的那些经典音频可视化器中不起作用。我也尝试了 matplotlib 的 pylab 的频谱图，但结果是一样的。
import os
import wave
import pylab
import math
from numpy import amax
from numpy import amin
def get_wav_info(wav_file,mi,mx):
wav = wave.open(wav_file, 'r')
frames = wav.readframes(-1)
sound_info = pylab.fromstring(frames, 'Int16')
frame_rate = wav.getframerate()
wav.close()
spectrum, freqs, t, im = pylab.specgram(sound_info, NFFT=1024, Fs=frame_rate)
n = 0
while n < 20:
for index,power in enumerate(spectrum[n]):
print("%s,%s,%s" % (n,int(round(t[index]*1000)),math.ceil(power*100)/100))
n += 1
get_wav_info("wave.wav",1,20)
关于如何获得可用于可视化的 dB 的任何提示？
基本上，我们显然已经从上面的代码中得到了我们需要的一切，只是如何使它 return 正常值？忽略 mi 和 mx，因为它们只是调整数组中的值以适应 mi..mx 间隔 - 这将用于可视化用途。如果我是正确的，此代码中的 spectrum returns 数组包含来自 freqs 数组的每个频率的振幅，根据 t 数组准时出现，但是这个值是如何工作的 - 如果它 return 是这些奇怪的值，它真的是振幅吗？如果是，如何将它转换为 dBs 例如。
tl;dr 我需要像音乐播放器那样的可视化工具输出，但它不应该实时工作，我只想要数据，但值不适合 wav 文件。
Edit2: 我注意到还有一个问题。对于 90 秒的 wav，t 数组包含直到 175.x 的时间，考虑到 frame_rate 对于 wav 文件是正确的，这看起来很奇怪。所以现在我们有 2 个问题：spectrum 似乎没有 return 正确的值（如果我们得到正确的时间，也许它会适合）并且 t 似乎 return 恰好加倍wav 的时间。
已修复：案例完全解决。
import os
import pylab
import math
from numpy import amax
from numpy import amin
from scipy.io import wavfile
frame_rate, snd = wavfile.read(wav_file)
sound_info = snd[:,0]
spectrum, freqs, t, im = pylab.specgram(sound_info,NFFT=1024,Fs=frame_rate,noverlap=5,mode='magnitude')
Specgram 需要一点调整，我只用 scipy.io 库（而不是 wave 库）加载了一个通道。同样，在没有将模式设置为幅度的情况下，它 return 是 10log10 而不是 20log10，这就是它没有 return 正确值的原因。
我在管理 FFT 数据时遇到了一点问题。我一直在寻找许多如何进行 FFT 的示例，但我无法从其中任何一个中得到我想要的东西。我有一个 44kHz 采样率的随机波形文件，我想每 X 毫秒获得 N 次谐波的幅度，假设 100 毫秒就足够了。我试过这段代码：
import scipy.io.wavfile as wavfile
import numpy as np
import pylab as pl
rate, data = wavfile.read("sound.wav")
t = np.arange(len(data[:,0]))*1.0/rate
p = 20*np.log10(np.abs(np.fft.rfft(data[:2048, 0])))
f = np.linspace(0, rate/2.0, len(p))
pl.plot(f, p)
pl.xlabel("Frequency(Hz)")
pl.ylabel("Power(dB)")
pl.show()
这是我使用的最后一个示例，我在 Whosebug 的某个地方找到了它。问题是，这得到了我想要的幅度，得到了频率，但根本没有时间。据我所知，FFT 分析是 3D 的，这是所有谐波的 "merged" 结果。我明白了：
X-axis = Frequency, Y-axis = Magnitude, Z-axis = Time (invisible)
根据我对代码的理解，t 是时间 - 看起来是这样，但代码中不需要 - 不过我们可能会需要它。 p 是幂（或幅度）数组，但它似乎是每个频率 f 的所有幅度的某个平均值，它是频率数组。我不想要 average/merged 值，我想要每 X 毫秒 N 次谐波的幅度。
长话短说，我们可以得到：所有频率的1级。
我们想要：N 个频率的所有幅度，包括出现特定幅度时的时间。
结果应类似于此数组：[时间、频率、振幅] 所以最后如果我们想要 3 个谐波，它看起来像：
[0,100,2.85489] #100Hz harmonic has 2.85489 amplitude on 0ms
[0,200,1.15695] #200Hz ...
[0,300,3.12215]
[100,100,1.22248] #100Hz harmonic has 1.22248 amplitude on 100ms
[100,200,1.58758]
[100,300,2.57578]
[200,100,5.16574]
[200,200,3.15267]
[200,300,0.89987]
不需要可视化，结果应该只是上面列出的数组（或 hashes/dictionaries）。
您似乎正在尝试实现 spectrogram，这是一系列功率谱估计，通常通过一系列（通常重叠的）FFT 来实现。由于您只有一个 FFT（频谱），因此您还没有时间维度。将 FFT 代码放入循环中，每次迭代处理一个样本块（例如 1024），连续块之间有 50% 的重叠。生成的频谱序列将是时间 v 频率 v 幅度的 3D 数组。
我不是 Python 人，但我可以给你一些伪代码，应该足以让你编码：
N = length of data input
N_FFT = no of samples per block (== FFT size, e.g. 1024)
i = 0 ;; i = index of spectrum within 3D output array
for block_start = 0 to N - block_start
block_end = block_start + N_FFT
get samples from block_start .. block_end
apply window function to block (e.g. Hamming)
apply FFT to windowed block
calculate magnitude spectrum (20 * log10( re*re + im*im ))
store spectrum in output array at index i
block_start += N_FFT / 2 ;; NB: 50% overlap
i++
end
根据@Paul R 的回答，scipy.signal.spectrogram 是 spectrogram function in scipy's signal processing module。
上面link的例子如下：
from scipy import signal
import matplotlib.pyplot as plt
# Generate a test signal, a 2 Vrms sine wave whose frequency linearly
# changes with time from 1kHz to 2kHz, corrupted by 0.001 V**2/Hz of
# white noise sampled at 10 kHz.
fs = 10e3
N = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
freq = np.linspace(1e3, 2e3, N)
x = amp * np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
#Compute and plot the spectrogram.
f, t, Sxx = signal.spectrogram(x, fs)
plt.pcolormesh(t, f, Sxx)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()
编辑： 哦，看起来这个 returns 值，但它们根本不适合音频文件。尽管它们可以用作频谱图上的幅度，但它们在您可以在许多音乐播放器中看到的那些经典音频可视化器中不起作用。我也尝试了 matplotlib 的 pylab 的频谱图，但结果是一样的。
import os
import wave
import pylab
import math
from numpy import amax
from numpy import amin
def get_wav_info(wav_file,mi,mx):
wav = wave.open(wav_file, 'r')
frames = wav.readframes(-1)
sound_info = pylab.fromstring(frames, 'Int16')
frame_rate = wav.getframerate()
wav.close()
spectrum, freqs, t, im = pylab.specgram(sound_info, NFFT=1024, Fs=frame_rate)
n = 0
while n < 20:
for index,power in enumerate(spectrum[n]):
print("%s,%s,%s" % (n,int(round(t[index]*1000)),math.ceil(power*100)/100))
n += 1
get_wav_info("wave.wav",1,20)
关于如何获得可用于可视化的 dB 的任何提示？
基本上，我们显然已经从上面的代码中得到了我们需要的一切，只是如何使它 return 正常值？忽略 mi 和 mx，因为它们只是调整数组中的值以适应 mi..mx 间隔 - 这将用于可视化用途。如果我是正确的，此代码中的 spectrum returns 数组包含来自 freqs 数组的每个频率的振幅，根据 t 数组准时出现，但是这个值是如何工作的 - 如果它 return 是这些奇怪的值，它真的是振幅吗？如果是，如何将它转换为 dBs 例如。
tl;dr 我需要像音乐播放器那样的可视化工具输出，但它不应该实时工作，我只想要数据，但值不适合 wav 文件。
Edit2: 我注意到还有一个问题。对于 90 秒的 wav，t 数组包含直到 175.x 的时间，考虑到 frame_rate 对于 wav 文件是正确的，这看起来很奇怪。所以现在我们有 2 个问题：spectrum 似乎没有 return 正确的值（如果我们得到正确的时间，也许它会适合）并且 t 似乎 return 恰好加倍wav 的时间。
已修复：案例完全解决。
import os
import pylab
import math
from numpy import amax
from numpy import amin
from scipy.io import wavfile
frame_rate, snd = wavfile.read(wav_file)
sound_info = snd[:,0]
spectrum, freqs, t, im = pylab.specgram(sound_info,NFFT=1024,Fs=frame_rate,noverlap=5,mode='magnitude')
Specgram 需要一点调整，我只用 scipy.io 库（而不是 wave 库）加载了一个通道。同样，在没有将模式设置为幅度的情况下，它 return 是 10log10 而不是 20log10，这就是它没有 return 正确值的原因。