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

这是一个非对数尺度的图，展示了一个合成波形文件，其中包含500Hz（逐渐减弱）和200Hz的正弦波，是用Goldwave制作的。

Answer 1

我对你的问题了解得不够，没法给出具体的答案。

不过根据我写快速傅里叶变换（FFT）的经验，这里有几个建议可以试试：

• 确保你遵循奈奎斯特定律。
• 如果你在查看FFT的线性输出，可能会很难看到自己的信号，觉得一切都坏了。确保你查看的是FFT幅度的分贝（dB）。比如可以用这个命令：“plot(10*log10(abs(fft(x))))”。
• 为你的FFT()函数创建一个单元测试，使用生成的数据，比如一个纯音。然后把同样的数据输入到Matlab的FFT()中。对比这两个输出数据的绝对值差异，确保最大绝对值差异大约是10^-6（也就是说，唯一的差异是由于小的浮点误差造成的）。
• 确保你在对数据进行窗函数处理。

如果这三件事都没问题，那么你的FFT应该是正常的，问题可能出在输入数据上。

• 检查输入数据，看看是否有削波现象 http://www.users.globalnet.co.uk/~bunce/clip.gif。

时间域的削波现象在频域中会表现为信号的镜像，且在特定的规律间隔处幅度较小。

Answer 2

快速傅里叶变换（FFT）因为是经过窗口处理和采样的，所以会导致别名现象（aliasing），而且在频率域中也会出现采样问题。简单来说，在时间域中进行滤波其实就是在频率域中进行乘法运算，所以你可以直接应用一个滤波器，这个滤波器就是把每个频率乘以一个特定的值。例如，在通过频带内的频率乘以1，而在其他地方乘以0。那些意外出现的值很可能是因为别名现象造成的，也就是高频信号被“折叠”到你看到的低频信号上。为了避免别名现象，原始信号的频率范围必须限制在你采样率的一半以内，否则就会出现别名现象。更需要关注的是，别名现象可能会扭曲你关心的频率范围，因为在这个频带内，你希望知道的频率是你预期的那个。

另外要记住的是，当你从一个波形文件中提取数据时，实际上是在进行数学上的平方波乘法。这会导致一个sinx/x的函数与频率响应进行卷积。为了减少这种影响，你可以用类似汉宁窗（Hanning window）的东西去乘以原始的窗口信号。

Answer 3

模拟的波形不应该像你的图那样显示FFT，所以肯定有什么地方出错了，问题可能不在FFT上，而是在输入的波形上。你图中的主要问题不是波动，而是1000赫兹附近的谐波和500赫兹的亚谐波。模拟波形不应该出现这些（比如，看看我下面的图）。

首先，你可能想先画出原始波形，这样可能会很明显地指向问题所在。而且，把波形解包成无符号短整型（即"H"）似乎有点奇怪，尤其是之后没有一个大的零频率成分。

我通过对波形进行裁剪，得到了一个与您的FFT非常接近的结果，这个裁剪是由亚谐波和更高的谐波（还有Trevor）建议的。你可能在模拟或解包时引入了裁剪。无论如何，我通过在numpy中创建波形来避免这个问题。

这是正确的FFT应该是什么样子的（基本上是完美的，除了由于窗函数导致的峰值扩展）

这是一个被裁剪的波形（与您的FFT非常相似，从亚谐波到1000赫兹附近的三个高谐波的精确模式）

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

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