filtering - DSP - 通过 FFT 在频域中进行滤波

Question

我一直在尝试 FFT 的外皮质实现，但遇到了一些问题。

每当我在调用 iFFT 之前修改频率箱的幅度时，生成的信号都会包含一些咔嗒声和爆裂声，特别是当信号中存在低频时(如鼓或贝司)。但是，如果我以相同的因子衰减所有箱，则不会发生这种情况。

让我举一个 4 样本 FFT 的输出缓冲区的示例:

// Bin 0 (DC)
FFTOut[0] = 0.0000610351563
FFTOut[1] = 0.0

// Bin 1
FFTOut[2] = 0.000331878662
FFTOut[3] = 0.000629425049

// Bin 2
FFTOut[4] = -0.0000381469727
FFTOut[5] =  0.0

// Bin 3, this is the first and only negative frequency bin.
FFTOut[6] =  0.000331878662
FFTOut[7] = -0.000629425049

输出由成对的 float 组成，每个 float 代表单个 bin 的实部和虚部。因此，bin 0(数组索引 0、1)将表示直流频率的实部和虚部。正如你所看到的，bin 1 和 3 都有相同的值(除了 Im 部分的符号)，所以我猜 bin 3 是第一个负频率，最后索引 (4, 5) 将是最后一个正频率频率仓。

然后为了衰减频率 bin 1，我会这样做:

// Attenuate the 'positive' bin
FFTOut[2] *= 0.5;
FFTOut[3] *= 0.5;

// Attenuate its corresponding negative bin.
FFTOut[6] *= 0.5;
FFTOut[7] *= 0.5;

对于实际测试，我使用 1024 长度的 FFT，并且我始终提供所有样本，因此不需要 0 填充。

// Attenuate
var halfSize = fftWindowLength / 2;
float leftFreq = 0f;
float rightFreq = 22050f; 
for( var c = 1; c < halfSize; c++ )
{
    var freq = c * (44100d / halfSize);

    // Calc. positive and negative frequency indexes.
    var k = c * 2;
    var nk = (fftWindowLength - c) * 2;

    // This kind of attenuation corresponds to a high-pass filter.
    // The attenuation at the transition band is linearly applied, could
    // this be the cause of the distortion of low frequencies?
    var attn = (freq < leftFreq) ? 
                    0 : 
                    (freq < rightFreq) ? 
                        ((freq - leftFreq) / (rightFreq - leftFreq)) :
                        1;

    // Attenuate positive and negative bins.
    mFFTOut[ k ] *= (float)attn;
    mFFTOut[ k + 1 ] *= (float)attn;
    mFFTOut[ nk ] *= (float)attn;
    mFFTOut[ nk + 1 ] *= (float)attn;
}

显然我做错了什么，但不知道是什么。

我不想使用 FFT 输出作为生成一组 FIR 系数的方法，因为我正在尝试实现一个非常基本的动态均衡器。

频域滤波的正确方法是什么？我缺少什么？

此外，是否真的需要衰减负频率？我见过 FFT 实现，但结果是否定的。频率值在合成前归零。

提前致谢。

Answer 1

有两个问题:使用 FFT 的方式以及特定的滤波器。

传统上，滤波是通过时域卷积来实现的。您是对的，将输入信号和滤波器信号的频谱相乘是等效的。但是，当您使用离散傅里叶变换 (DFT)(通过快速傅里叶变换算法实现以提高速度)时，您实际上计算的是真实频谱的采样版本。这有很多含义，但与滤波最相关的含义是时域信号是周期性的。

这是一个例子。考虑一个周期为 1.5 个周期的正弦输入信号 x，以及一个简单的低通滤波器 h。在 Matlab/Octave 语法中:

N = 1024;
n = (1:N)'-1; %'# define the time index
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4
h = [h./sum(h)]; %# normalize DC gain

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra
y = real(y); %# due to numerical error, y has a tiny imaginary part
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here
plot(y);

这是图表:

block 开头的故障根本不是我们所期望的。但如果你考虑一下fft(x)，它就有意义了。离散傅里叶变换假设信号在变换 block 内是周期性的。据 DFT 所知，我们要求对其进行一个周期的变换:

这导致了使用 DFT 进行过滤时的第一个重要考虑因素:您实际上正在实现 circular convolution ，不是线性卷积。因此，当您考虑数学时，第一张图中的“故障”并不是真正的故障。那么问题就变成了:有没有办法解决周期性问题？答案是肯定的:使用overlap-save processing 。本质上，您可以像上面那样计算 N 长乘积，但只保留 N/2 个点。

Nproc = 512;
xproc = zeros(2*Nproc,1); %# initialize temp buffer
idx = 1:Nproc; %# initialize half-buffer index
ycorrect = zeros(2*Nproc,1); %# initialize destination
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer
    idx = idx + Nproc; %# step half-buffer index
end

这是y Correct的图表:

这张图很有意义 - 我们预计滤波器会出现启动 transient ，然后结果稳定为稳态正弦响应。请注意，现在 x 可以是任意长。限制为 Nproc > 2*min(length(x),length(h))。

现在讨论第二个问题:特定的过滤器。在循环中，您创建一个滤波器，其频谱本质上是 H = [0 (1:511)/512 1 (511:-1:1)/512]'; 如果您这样做 hraw = real(ifft(H));绘制(hraw)，你得到:

很难看出，但在图表的最左边缘有一堆非零点，然后在最右边缘还有一堆非零点。使用 Octave 的内置 freqz 函数查看我们看到的频率响应(通过执行 freqz(hraw) ):

幅度响应有很多从高通包络到零的纹波。同样，DFT 固有的周期性也在发挥作用。就 DFT 而言，hraw 一遍又一遍地重复。但是，如果您采用一个周期的 hraw，就像 freqz 那样，它的频谱与周期版本的频谱有很大不同。

所以让我们定义一个新信号:hrot = [hraw(513:end) ; hraw(1:512)]; 我们只需旋转原始 DFT 输出以使其在 block 内连续。现在让我们使用 freqz(hrot) 看一下频率响应:

好多了。所需的信封就在那里，没有任何波纹。当然，现在的实现并不那么简单，你必须用 fft(hrot) 进行完整的复数乘法，而不是仅仅缩放每个复数箱，但至少你会得到正确的答案.

请注意，为了速度，您通常会预先计算填充的 h 的 DFT，我将其单独留在循环中，以便更轻松地与原始值进行比较。