削波 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

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

---

编辑：这里是基数的示例-2 FFT 在纯 Python 中实现，通过一个简单的测试例程来查找矩形脉冲的 FFT，[1.,1.,1.,1.,0.,0.,0.,0.]< /代码>。 您可以在 codepad 上运行该示例，并看到该序列的 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，然后到最高正频率。

Answer 2

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

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

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

Answer 4

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

但根据我自己编写 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

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