裁剪 FFT 矩阵

Question

音频处理对我来说相当新。目前使用 Python Numpy 处理波形文件。在计算 FFT 矩阵后，我得到了不存在频率的噪声功率值。我对可视化数据感兴趣，准确性不是一个高优先级。是否有一种安全的方法来计算裁剪值以删除这些值，或者我应该使用每个样本集的所有 FFT 矩阵来得出一个平均数？

问候

编辑：

    from numpy import *
    import wave
    import pymedia.audio.sound as sound
    import time, struct
    from pylab import ion, plot, draw, show

    fp = wave.open("500-200f.wav", "rb")
    sample_rate = fp.getframerate()
    total_num_samps = fp.getnframes()
    fft_length = 2048.
    num_fft = (total_num_samps / fft_length ) - 2
    temp = zeros((num_fft,fft_length), float)

    for i in range(num_fft):
        tempb = fp.readframes(fft_length);
        data = struct.unpack("%dH"%(fft_length), tempb)
        temp[i,:] = array(data, short) 
    pts = fft_length/2+1
    data = (abs(fft.rfft(temp, fft_length)) / (pts))[:pts]

    x_axis = arange(pts)*sample_rate*.5/pts
    spec_range = pts
    plot(x_axis, data[0])
    show()

这是使用 Goldwave 创建的包含 500hz（淡出）+ 200hz 正弦波的合成波文件的非对数比例图。

Answer 1

模拟波形不应该像您的图那样显示 FFT，所以有些地方很不对劲，可能不是 FFT，而是输入波形。您绘图中的主要问题不是波纹，而是 1000 Hz 左右的谐波和 500 Hz 的次谐波。模拟波形不应显示任何这些（例如，请参见下面的图表）。

首先，您可能只想尝试绘制原始波形，这可能会指出一个明显的问题。此外，将波解包到无符号短裤（即“H”）似乎很奇怪，尤其是在此之后没有大的零频率分量。

正如分谐波和高次谐波（和 Trevor）所建议的那样，通过对波形应用削波，我能够得到与您的 FFT 非常接近的副本。您可以在模拟或解包中引入剪裁。无论哪种方式，我通过在 numpy 中创建波形来绕过这一点。

这就是正确的 FFT 应该是什么样子（即基本上是完美的，除了由于开窗而使峰变宽）

这是一个被削波的波形（与您的 FFT 非常相似，从次谐波到 1000 Hz 左右三个高次谐波的精确模式）

这是我用来生成这些的代码

from numpy import *
from pylab import ion, plot, draw, show, xlabel, ylabel, figure

sample_rate = 20000.
times = arange(0, 10., 1./sample_rate)
wfm0 = sin(2*pi*200.*times)
wfm1 = sin(2*pi*500.*times) *(10.-times)/10.
wfm = wfm0+wfm1
#  int test
#wfm *= 2**8
#wfm = wfm.astype(int16)
#wfm = wfm.astype(float)
#  abs test
#wfm = abs(wfm)
#  clip test
#wfm = clip(wfm,  -1.2, 1.2)

fft_length = 5*2048.
total_num_samps = len(times)
num_fft = (total_num_samps / fft_length ) - 2
temp = zeros((num_fft,fft_length), float)

for i in range(num_fft):
    temp[i,:] = wfm[i*fft_length:(i+1)*fft_length] 
pts = fft_length/2+1
data = (abs(fft.rfft(temp, fft_length)) / (pts))[:pts]

x_axis = arange(pts)*sample_rate*.5/pts
spec_range = pts
plot(x_axis, data[2], linewidth=3)
xlabel("freq (Hz)")
ylabel('abs(FFT)')
show()

Answer 2

FFT 因为它们是加窗的，sampled 导致 aliasing 和频域采样。时域中的过滤只是频域中的乘法，因此您可能只想应用一个过滤器，该过滤器只是将每个频率乘以您正在使用的过滤器的函数值。例如，在通带中乘以 1，在其他区域中乘以零。意外的值可能是由混叠引起的，其中较高的频率被折叠到您所看到的频率。 original signal needs to be band limited to half your sampling rate 或者你会得到aliasing。更令人担忧的是扭曲感兴趣区域的混叠，因为对于这个频带，您想知道频率来自预期的频率。

要记住的另一件事是，当您从波形文件中获取一条数据时，您在数学上将其乘以方波。这会导致 sinx/x 与频率响应进行卷积以将其最小化，您可以将原始加窗信号与 Hanning window 之类的信号相乘。

Answer 3

对于一维 FFT，值得一提的是，第一个元素（索引 [0]）包含 DC（零频率）项，元素 [1:N/2] 包含正频率，元素 [N/2+1:N-1] 包含负频率。由于您没有提供有关 FFT 输出的代码示例或其他信息，因此我不能排除“不存在频率下的噪声功率值”不仅仅是频谱的负频率的可能性。

---

编辑：Here 是一个用纯 Python 实现的 radix-2 FFT 示例，带有一个简单的测试例程，可以找到矩形脉冲的 FFT，[1.,1.,1.,1.,0.,0.,0.,0.]。您可以在codepad 上运行示例，并看到该序列的 FFT 是

[0j,                    Negative frequencies
(1+0.414213562373j),    ^
0j,                     |
(1+2.41421356237j),     |
(4+0j),                <= DC term
(1-2.41421356237j),     |
0j,                     v
(1-0.414213562373j)]    Positive frequencies

请注意，代码按频率升序打印傅立叶系数，即从最高负频率到 DC，然后再到最高正频率。

Answer 4

我对你的问题了解得不够多，无法实际回答任何具体问题。

但根据我自己编写 FFT 的经验，可以尝试以下几点：

• 确保您遵循奈奎斯特规则
• 如果您正在查看 FFT 的线性输出...您将无法查看自己的信号并认为一切都已损坏。确保您正在查看 FFT 幅度的 dB。 （即“情节（10 * log10（abs（fft（x））））”）
• 通过像纯音一样输入生成的数据，为您的 FFT() 函数创建一个 unitTest。然后将相同的生成数据提供给 Matlab 的 FFT()。在两个输出数据系列之间进行绝对值差异，并确保最大绝对值差异类似于 10^-6（即唯一的差异是由小浮点错误引起的）
• 确保您是windowing your data

如果这三件事都有效，那么您的 fft 就可以了。而您的输入数据可能就是问题所在。

• 查看输入数据是否有clipping http://www.users.globalnet.co.uk/~bunce/clip.gif

时域削波显示为频域中信号的镜像，以特定的规则间隔呈现，幅度较小。