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

Question

我一直在玩 FFT 的 Exocortex 实现，但遇到了一些问题。

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

让我举一个 4-sample 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

输出由成对的浮点数组成，每个浮点数代表单个 bin 的实部和虚部。因此，bin 0（数组索引 0、1）将代表 DC 频率的实部和虚部。如您所见，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 实现，其中 neg。频率值在合成前归零。

提前致谢。

Answer 1

首先，关于标准化：这是一个已知（非）问题。 DFT/IDFT 需要一个因子 1/sqrt(N)（除了标准的余弦/正弦因子）（正反），以使它们对称且真正可逆。另一种可能性是将其中一个（正向或反向）除以 N，这是方便和品味的问题。 FFT 例程通常不执行这种归一化，用户应该知道它并按照他的喜好进行归一化。 See

第二：在（比如说）16 点 DFT 中，您所说的 bin 0 将对应于零频率 (DC)，bin 1 低频.. . bin 4 中频，bin 8 到最高频率，bins 9...15 到“负频率”。那么，在您的示例中， bin 1 实际上既是低频又是中频。除了这个考虑之外，您的“均衡”在概念上没有任何错误。我不明白您所说的“信号在低频下失真”是什么意思。你是怎么观察到的？

Answer 2

您是否将直流频率样本的值衰减到零？在您的示例中，您似乎根本没有减弱它。由于您正在实现高通滤波器，因此您还需要将 DC 值设置为零。

这可以解释低频失真。如果 DC 值不为零，因为过渡较大，低频时的频率响应中会出现很多纹波。

以下是 MATLAB/Octave 中的示例，用于演示可能发生的情况：

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,4) linspace(1,0,8) zeros(1,3) 1 zeros(1,4) linspace(0,1,8) ones(1,4)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude');

请注意，在我的代码中，我创建了一个 DC 值不为零的示例，然后是突然变为零，然后是斜升。然后我将IFFT转换为时域。然后我对该时域信号执行零填充 fft（当您传入大于输入信号的 fft 大小时，由 MATLAB 自动完成）。时域中的零填充导致频域中的插值。使用它，我们可以看到过滤器将如何在过滤器样本之间做出响应。

要记住的最重要的事情之一是，即使您通过衰减 DFT 的输出来设置给定频率的滤波器响应值，这也不能保证采样点之间出现的频率。这意味着您的更改越突然，样本之间的过冲和振荡就会越多。

现在回答您关于如何进行此过滤的问题。有很多方法，但最容易实现和理解的方法之一是窗口设计方法。您当前设计的问题是过渡宽度很大。大多数时候，你会希望尽可能快的过渡，尽可能少的波动。

在接下来的代码中，我将创建一个理想的过滤器并显示响应：

N = 32;
os = 4;
Fs = 1000;
X = [ones(1,8) zeros(1,16) ones(1,8)];
x = ifftshift(ifft(X));
Xos = fft(x, N*os);
f1 = linspace(-Fs/2, Fs/2-Fs/N, N);
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os);

hold off;
plot(f2, abs(Xos), '-o');
hold on;
grid on;
plot(f1, abs(X), '-ro');
hold off;
xlabel('Frequency (Hz)');
ylabel('Magnitude'); 

请注意，突然变化会引起很多振荡。

FFT 或离散傅里叶变换是傅里叶变换的采样版本。傅里叶变换应用于连续范围内的信号——无穷大到无穷大，而 DFT 应用于有限数量的样本。当使用 DFT 时，这实际上导致了时域中的平方窗口（截断），因为我们只处理有限数量的样本。不幸的是，方波的 DFT 是一个 sinc 类型的函数 (sin(x)/x)。

过滤器中的急剧转换（一个样本中从 0 到 1 的快速跳转）的问题在于，它在时域中具有非常长的响应，该响应被方形窗口截断。因此，为了帮助最小化这个问题，我们可以将时域信号乘以一个更渐变的窗口。如果我们通过添加以下行来乘以汉宁窗：

x = x .* hanning(1,N).';

在接受 IFFT 后，我们得到以下回复：

所以我建议尝试实现窗口设计方法，因为它相当简单（有更好的方法，但它们变得更复杂）。由于您正在实现均衡器，我假设您希望能够动态更改衰减，因此我建议在参数发生变化时计算并将滤波器存储在频域中，然后您可以应用它通过获取输入缓冲区的 fft，乘以频域滤波器样本，然后执行 ifft 以返回时域，到每个输入音频缓冲区。这将比您为每个样本执行的所有分支更有效。

Answer 3

有两个问题：您使用 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);

这是图表：

块开头的故障根本不是我们所期望的。但是，如果您考虑fft(x)，这是有道理的。离散傅里叶变换假设信号在变换块内是周期性的。据 DFT 所知，我们要求转换其中一个周期：

这导致了使用 DFT 进行过滤时的第一个重要考虑因素：您实际上是在实现 circular convolution，而不是线性卷积。因此，当您考虑数学时，第一张图中的“故障”并不是真正的故障。那么问题就变成了：有没有办法解决周期性问题？答案是肯定的：使用overlap-save processing。本质上，您如上所述计算 N-long 产品，但只保留 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

这是ycorrect 的图表：

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

现在进入第二个问题：特定过滤器。在您的循环中，您创建一个过滤器，其频谱本质上为H = [0 (1:511)/512 1 (511:-1:1)/512]'; 如果您执行hraw = real(ifft(H)); plot(hraw)，您将获得：

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

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

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

好多了。所需的信封就在那里，没有所有的涟漪。当然，现在的实现并不是那么简单，你必须做一个完整的复数乘以fft(hrot)，而不是仅仅缩放每个复数 bin，但至少你会得到正确的答案。

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

Answer 4

您的主要问题是在短时间内没有很好地定义频率。对于低频尤其如此，这就是为什么您最注意到那里的问题的原因。

因此，当您从声音序列中取出非常短的片段，然后对它们进行过滤时，过滤后的片段不会以产生连续波形的方式进行过滤，您会听到片段之间的跳跃，这就是产生点击这里。

例如，取一些合理的数字：我从 27.5 Hz 的波形（钢琴上的 A0）开始，以 44100 Hz 数字化，它看起来像这样（红色部分是 1024 个样本长）：

所以首先我们将从 40Hz 的低通开始。因此，由于原始频率小于 40Hz，具有 40Hz 截止频率的低通滤波器应该不会产生任何影响，我们将获得几乎与输入完全匹配的输出。正确的？ 错了，错了，错了——这基本上是你问题的核心。问题在于，对于较短的部分，27.5 Hz 的idea 没有明确定义，并且不能在 DFT 中很好地表示。

通过查看下图中的 DFT，可以看出 27.5 Hz 在短段中并不是特别有意义。请注意，虽然较长段的 DFT（黑点）在 27.5 Hz 处显示出峰值，但较短段（红点）则没有。

显然，然后在 40Hz 以下进行滤波，只会捕获 DC 偏移，而 40Hz 低通滤波器的结果如下图绿色所示。

蓝色曲线（采用 200 Hz 截止频率）开始匹配得更好。但请注意，使它匹配得很好的不是低频，而是高频的包含。直到我们在短段中包含所有可能的频率，最高可达 22KHz，我们才能最终获得原始正弦波的良好表示。

所有这一切的原因是 27.5 Hz 正弦波的一小部分不是 27.5 Hz 正弦波，它的 DFT 与 27.5 Hz 没有太大关系。