裁剪FFT矩阵

Question

音频处理对我来说是个很新的东西。目前正在使用Python Numpy处理wave文件。在计算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创建的包含500 is (淡出)+200 is正弦波的合成波文件的非对数比例图。

​

​

Answer 1

模拟的波形不应该显示像您的图那样的FFT，所以有些东西是非常错误的，可能不是FFT，而是输入波形。你图中的主要问题不是波纹，而是1000赫兹附近的谐波和500赫兹的次谐波。模拟的波形不应该显示任何这些(例如，请参见下面的图表)。

首先，您可能只想尝试绘制出原始波形，这可能会指出一个明显的问题。此外，有一个波解包到未签名的短裤似乎很奇怪，即。"H"，尤其是在此之后不会有很大的零频分量。

通过对波形应用裁剪，我能够得到与你的FFT非常接近的副本，正如亚谐波和高次谐波(以及Trevor)所建议的那样。您可以在模拟或解包过程中引入剪裁。无论哪种方式，我都通过在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

快速傅立叶变换是加窗的，sampled也会在频域中产生aliasing和采样。时域中的滤波只是频域中的乘法，因此您可能只想应用一个滤波器，它只是将每个频率乘以您正在使用的滤波器的函数的值。例如，在通带中乘以1，在其他情况下乘零。意外的值可能是由混叠引起的，其中较高的频率被折叠到您所看到的频率。original signal needs to be band limited to half your sampling rate，否则您将获得aliasing。更值得关注的是混叠，它扭曲了感兴趣的区域，因为对于这段频率，您希望知道该频率来自于预期的频率。

另一件要记住的事情是，当你从wave文件中抓取一段数据时，你在数学上将它乘以一个方波。这会导致sinx/x与频率响应卷积，为了最小化这一点，您可以将原始加窗信号与Hanning window之类的东西相乘。

Answer 3

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

Python ：Here是一个用纯实现的基数-2FFT的例子，它有一个简单的测试例程，可以找到矩形脉冲[1.,1.,1.,1.,0.,0.,0.,0.]的FFT。您可以在codepad上运行该示例，并看到该序列的快速傅立叶变换为

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

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