【问题标题】:any rules of thumb how to smooth FFT spectrum to prevent artifacts when hand-tweaking?任何经验法则如何在手动调整时平滑 FFT 频谱以防止伪影?
【发布时间】:2011-04-19 03:31:12
【问题描述】:

我有一个 FFT 幅度谱,我想从中创建一个滤波器,选择性地通过周期性噪声源(例如正弦波杂散)并将与随机背景噪声相关的频率区间归零。我知道一旦这个过滤器IFFT回到时域,频率域中的急剧转换会产生振铃伪影......所以我想知道是否有任何经验法则如何平滑这种过滤器中的转换以避免这种情况铃声。

例如,如果 FFT 有 1M 频率区间,并且有五个杂散从背景噪声基底中伸出,我想将除与五个杂散中的每一个相关联的峰值区间之外的所有区间归零。问题是如何处理相邻的支线箱以防止时域中的伪影。例如,是否应将支线箱每一侧的箱设置为 50% 幅度?是否应该在支线箱的两侧使用两个箱(最近的一个在 50%,下一个在 25% 等)?任何想法都非常感谢。谢谢!

【问题讨论】:

    标签: fft frequency dft


    【解决方案1】:

    我喜欢下面的方法:

    • 创建理想的幅度谱(记住使其关于 DC 对称)
    • 时域逆变换
    • 将块旋转一半的块大小
    • 应用 Hann 窗口

    我发现它可以创建相当平滑的频域结果,尽管我从未尝试过像您建议的那样尖锐的东西。您可能可以通过使用 Kaiser-Bessel 窗口来制作更清晰的过滤器,但您必须适当地选择参数。通过更锐利,我猜也许你可以将旁瓣减少 6 dB 左右。

    这里是一些示例 Matlab/Octave 代码。为了测试结果,我使用了freqz(h, 1, length(h)*10);

    function [ht, htrot, htwin] = ArbBandPass(N, freqs)
    %# N = desired filter length
    %# freqs = array of frequencies, normalized by pi, to turn into passbands
    %# returns raw, rotated, and rotated+windowed coeffs in time domain
    
    if any(freqs >= 1) || any(freqs <= 0)
        error('0 < passband frequency < 1.0 required to fit within (DC,pi)')
    end
    
    hf = zeros(N,1); %# magnitude spectrum from DC to 2*pi is intialized to 0
    %# In Matlabs FFT, idx 1 -> DC, idx 2 -> bin 1, idx N/2 -> Fs/2 - 1, idx N/2 + 1 -> Fs/2, idx N -> bin -1
    idxs = round(freqs * N/2)+1; %# indeces of passband freqs between DC and pi
    hf(idxs) = 1; %# set desired positive frequencies to 1
    hf(N - (idxs-2)) = 1; %# make sure 2-sided spectrum is symmetric, guarantees real filter coeffs in time domain
    ht = ifft(hf); %# this will have a small imaginary part due to numerical error
    if any(abs(imag(ht)) > 2*eps(max(abs(real(ht)))))
        warning('Imaginary part of time domain signal surprisingly large - is the spectrum symmetric?')
    end
    ht = real(ht); %# discard tiny imag part from numerical error
    htrot = [ht((N/2 + 1):end) ; ht(1:(N/2))]; %# circularly rotate time domain block by N/2 points
    win = hann(N, 'periodic'); %# might want to use a window with a flatter mainlobe
    htwin = htrot .* win;
    htwin = htwin .* (N/sum(win)); %# normalize peak amplitude by compensating for width of window lineshape
    

    【讨论】:

    • 非常感谢 mtrw。我正在考虑手动调整相邻 bin 上的 FFT 系数,但似乎您正在使用 Hann 窗口来实现类似的效果(正确吗?伪影将出现在时域中,而 Hann 窗口会以某种方式改善它) .在我的情况下,“理想幅度谱”将是二进制表示,其中(继续上面的示例),我有五个 1 和其余的零 - 这也是你在这里推荐的吗?另外,我不确定“将块旋转一半块大小”是什么意思。什么是块?能举个例子吗?
    • 另外,以 DC 为中心的频谱有什么优势?
    • 快速回答: 1. 是的,我建议窗口负责处理相邻的垃圾箱。 2. 旋转块意味着在逆 FFT 之后交换时域滤波器的前半部分和后半部分。 Blocksize 是过滤器中的点数。 3. 您需要确保逆 FFT 提供实时时域系数。光谱必须是 Hermitian 对称的(第二个 N/2 点是第一个 N/2 的复共轭)。检查您的 FFT 文档以了解他们希望您如何安排这些点。有时间我会创建一个示例...
    • 谢谢mtrw。我有一个包含 N 个真实数据点的波形。我正在从前 1...N/2 个点构建一个过滤器,因此过滤器长度为 M=N/2。使用 Matlab 的 fft 例程。将过滤器和数据数组都归零填充到 M+N-1 以允许消耗卷积空间。只对稳态瞬变感兴趣,并考虑使用窗口来帮助马刺脱颖而出,但仍在争论中。
    • 我添加了一个代码示例,该示例生成一个 N 点时域滤波器,其频率响应具有任意数量的 bin 传递能量。但我有一些问题: 1. 我完全不明白你对“稳态瞬变”的评论。 2. 我也不明白你打算如何使用过滤器。通常,过滤器比应用它们的数据要短。您需要多少过渡带宽?
    猜你喜欢
    • 2020-08-30
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-07-24
    • 1970-01-01
    • 2019-08-05
    • 2022-12-20
    • 2021-09-25
    相关资源
    最近更新 更多