python - 削波 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()