DSP-通过FFT在频域中过滤

Question

我一直在使用FFT的Exocortex实现来玩耍，但是我遇到了一些问题。

每当我在调用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

输出由一对浮子组成，每个浮子代表单个垃圾箱的真实和想象中的部分。因此，bin 0（阵列索引0，1）将代表直流频率的真实和虚构部分。如您所见，垃圾箱1和3都具有相同的值（除了IM部分的符号外），所以我猜Bin 3是第一个负频率，最后索引（4，5）将是最后一个正阳性频率箱。

然后，要衰减频率箱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）（使用快速傅立叶变换算法实现）时，您实际上计算了真实频谱的采样版本。这有很多含义，但是与过滤最相关的是，时间域信号是周期性的。

这是一个例子。考虑正弦输入信号 x 在此期间有1.5个周期，并且简单的低通滤波器 h. 。在MATLAB/八度语法中：

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过滤时的第一个重要考虑：您实际上正在实施 圆形卷积, ，不是线性卷积。因此，当您考虑数学时，第一个图中的“小故障”并不是真正的故障。那么问题就变成了：有没有办法解决周期性？答案是肯定的：使用 重叠保存处理. 。从本质上讲，您按照上述计算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), ， 你得到：

很难看到，但是图表的最左边缘有很多非零点，然后在最右边的边缘有更多。使用八度的内置 freqz 查看我们看到的频率响应的功能（通过 freqz(hraw)):

幅度响应从高通信封下降到零有很多涟漪。同样，DFT固有的周期性正在起作用。就DFT而言 hraw 一遍又一遍地重复。但是，如果您花了一段时间 hraw, ， 作为 freqz 确实，它的频谱与定期版本完全不同。

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

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

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

Answer 2

您的主要问题是 频率在短时间间隔不正确. 。对于低频尤其如此，这就是为什么您注意到那里最问题的原因。

因此，当您从声音列车中取出真正的短段，然后将其过滤时，过滤后的段不会以产生连续波形的方式过滤，并且您会听到片段之间的跳跃，这就是在此处产生点击的原因。

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

因此，首先，我们将从40Hz的低通票开始。因此，由于原始频率小于40Hz，因此具有40Hz截止的低通滤波器应真正没有任何效果，我们将获得几乎与输入完全匹配的输出。正确的？ 错，错，错了 - 这基本上是您问题的核心。问题在于，对于简短的部分 主意 在27.5 Hz中没有明确定义，在DFT中不能很好地表示。

通过查看下图中的DFT，可以看到27.5 Hz在短段中并没有特别有意义。请注意，尽管较长的段的DFT（黑点）显示为27.5 Hz的峰值，但短一个（红点）没有。

显然，然后将其过滤在40Hz以下，将仅捕获直流偏移，而40Hz低通滤波器的结果在下面以绿色显示。

蓝色曲线（以200 Hz的截止为止）开始变得更好。但是请注意，使其匹配的不是低频，而是高频的包含。直到我们在短段中包括每个可能的频率，最多22kHz，我们最终才能很好地表示原始的正弦波。

所有这一切的原因是，一小部分27.5 Hz正弦波为 不是 27.5 Hz正弦波，DFT与27.5 Hz无关。

Answer 3

您是否正在将直流频率样本的值降低到零？看来您在示例中根本没有减弱它。由于您正在实现高通滤波器，因此您也需要将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');

请注意，在我的代码中，我正在创建一个非零的直流值的示例，然后突然更改为零，然后逐渐增加。然后，我将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）在于，这在时域中的响应很长，该响应被方形窗口截断。因此，为了帮助最大程度地减少此问题，我们可以将时间域信号乘以更逐渐的窗口。如果我们通过添加该行将Hanning窗口乘以：

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

服用ifft后，我们得到了此回应：

因此，我建议尝试实现窗口设计方法，因为它非常简单（有更好的方法，但它们变得更加复杂）。由于您正在实现均衡器，因此我假设您希望能够随时更改衰减，因此我建议在参数发生更改时计算并将过滤器存储在频域中，然后您只需应用它到每个输入音频缓冲区，通过取输入缓冲区的FFT，乘以您的频域滤波器样品，然后执行iFFT以返回时间域。这将比您为每个样本所做的所有分支要高得多。

Answer 4

首先，关于标准化：这是一个已知的（非）问题。 DFT/IDFT需要一个因素 1/sqrt（n） （除了标准的余弦/正弦因素外）在每个（直接逆）中，使它们具有相似和真正的可逆性。另一种可能性是将其中一个（直接或逆）除以 n, ，这是一个方便和品味的问题。通常，FFT例程不会执行此归一化，预计用户会意识到它并在他喜欢的情况下正常化。 看

第二：在（说）16点DFT中，您称之为 bin 0 将对应于零频率（DC）， bin 1 低频... 垃圾箱4 中等频率， 本8 达到最高频率和 垃圾箱9 ... 15 到“负频率”。在你的示例中， bin 1 实际上是低频和中频。除了这种考虑之外，您的“均等化”中没有任何概念上的错误。我不明白你的意思 “信号以低频扭曲”. 。您如何观察？