2

音频处理对我来说很新鲜。目前使用 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 正弦波的合成波文件的非对数刻度图。

4

4 回答 4

3

模拟波形不应该像你的图那样显示 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()
于 2009-06-01T01:27:00.210 回答
2

FFT 因为它们是加窗和采样的，所以也会导致频域中的混叠和采样。时域中的过滤只是频域中的乘法，因此您可能只想应用一个过滤器，该过滤器只是将每个频率乘以您正在使用的过滤器的函数值。例如，在通带中乘以 1，在其他区域中乘以零。意外的值可能是由混叠引起的，其中较高的频率被折叠到您所看到的频率。原始信号需要将频带限制为采样率的一半，否则会出现混叠. 更令人担忧的是扭曲感兴趣区域的混叠，因为对于这个频带，您想知道频率来自预期的频率。

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

于 2009-05-31T23:49:36.443 回答
1

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

编辑：这是一个用纯 Python 实现的 radix-2 FFT 示例，它有一个简单的测试例程，可以找到矩形脉冲的 FFT [1.,1.,1.,1.,0.,0.,0.,0.]，. 您可以在键盘上运行该示例并查看该序列的 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，然后再到最高正频率。

于 2009-06-01T04:24:43.613 回答
1

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

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

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

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

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

于 2009-06-01T13:34:16.277 回答