如何在Python中从FFT获取时间/频率

Question

我在处理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()

这是我使用的最后一个例子，在stackoverflow上找到的。问题是，它得到了我想要的大小和频率，但完全没有时间。据我所知，FFT分析是三维的，而这是所有谐波的“合并”结果。我得到了这个：X轴 = 频率、Y轴 = 大小、Z轴 = 时间（不可见）。

根据我的代码理解，t是时间，似乎是可以使用的，但在代码中不需要。p是功率数组（或大小），但看起来像是每个频率f的所有大小的平均值，而f是频率数组。我不想要平均/合并值，我想要每X毫秒N个谐波的幅度大小。

简而言之，我们可以得到：所有频率的1个幅度大小。

我们想要：N个频率的所有幅度大小，包括出现某个幅度大小的时间。

结果应该像这样的数组：[时间，频率，振幅]。因此，如果我们想要三个谐波，它会看起来像：

[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]

可视化并不是必需的，结果只需要按上面列出的数组（或哈希/字典）形式呈现。

- Dulcia

---

快速傅里叶变换（FFT）算法计算序列的离散傅里叶变换（DFT）或其反变换。傅里叶分析将信号从其原始域（通常是时间或空间）转换为频率域表示，反之亦然。我认为一旦在原始信号上应用了傅里叶变换，就不应该再考虑时间。它会被转换为频率域。同样地，当您在频率域信号上应用反傅里叶变换时，您会得到时间域信号。在此阅读更多信息。https://en.wikipedia.org/wiki/Fast_Fourier_transform - Sagar Waghmode

谢谢您的评论，虽然您向我解释了算法的工作原理，但我仍然不知道是否可能从中获得这样的输出，或者是否需要完全不同的方法。如果不用FFT，那么如何才能获得我描述的输出呢？了解一点FFT的工作原理并不能解决问题。 - Dulcia

嗯，如果我反变换傅里叶变换，我会得到时域信号，但那是原始信号，是吧？除此之外，我还不知道在哪里可以获得所有三个值。 - Dulcia

Answer 1

在@Paul R的回答之后，scipy.signal.spectrogram是scipy信号处理模块中的谱图函数。

上面链接中的示例如下：

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()

- Lee

---

从f、t和Sxx获取我需要的所有三个值非常容易。问题将是将wav文件导入其中，然后它应该可以顺利运行。然而，我认为scipy库将在所有组件中兼容，但似乎并非如此。scipy.io中的wavfile.read从wav文件中生成一个ndarray，但不能作为信号spectrogram的输入，即使上面的代码中的x也是ndarray。我完全没有头绪，因为文档似乎没有展示与scipy.io.wavfile.read的任何连接。 - Dulcia

Answer 2

看起来您正在尝试实现 频谱图(spectrogram)，它是一系列功率谱估计的序列，通常使用一系列（通常重叠的）FFT实现。由于您只有一个FFT（谱），因此尚未具备时间维度。将您的FFT代码放入循环中，并在每次迭代中处理一个样本块（例如1024），连续块之间具有50％的重叠。然后生成的频谱序列将是一个三维数组，包含时间、频率和幅度。

虽然我不是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

---

我知道你的意思，并且可以确认spectogram是我正在寻找的。然而，作为一个新手，我不知道如何做到这一点。有没有人能给我一些提示或完整的示例？ - Dulcia

Answer 3

编辑：哦，看起来这个返回值，但它们根本不适合音频文件。即使它们可以用作频谱图上的幅度，它们也不能在许多音乐播放器中看到的那些经典音频可视化器中使用。我还尝试了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值吗？ 基本上，我们显然从上面的代码中得到了所有需要的东西，只是如何使其返回正常值？忽略mi和mx，因为这些只是调整数组中的值以适应mi..mx间隔 - 这将用于可视化使用。如果我没错，此代码中的spectrum返回包含每个频率的幅度的数组的数组，这些频率根据t数组在时间上存在，但该值如何工作-它是否真的是幅度，如果它返回这些奇怪的值并且如果是，如何将其转换为例如dB。 简而言之，我需要像音乐播放器一样的可视化器输出，但它不应实时工作，我只想要数据，但值不适合wav文件。

编辑2：我注意到还有一个问题。对于90秒的wav文件，t数组包含的时间只到175.x，这似乎非常奇怪，考虑到frame_rate与wav文件正确匹配。所以现在我们有两个问题： spectrum似乎不能返回正确的值（也许如果我们得到正确的时间就可以解决），而t似乎返回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 库）。此外，如果没有将模式设置为幅度，它会返回 10log10 而不是 20log10，这就是为什么它没有返回正确值的原因。

- Dulcia