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 回答