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

发布于 2024-09-03 20:05:18 字数 2015 浏览 9 评论 0原文

我一直在尝试 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

输出由成对的浮点数组成，每个浮点数代表单个 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 实现，但结果是否定的。频率值在合成之前归零。

提前致谢。

原文

I've been playing around a little with the Exocortex implementation of the FFT, but I'm having some problems.

Whenever I modify the amplitudes of the frequency bins before calling the iFFT the resulting signal contains some clicks and pops, especially when low frequencies are present in the signal (like drums or basses). However, this does not happen if I attenuate all the bins by the same factor.

Let me put an example of the output buffer of a 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

The output is composed of pairs of floats, each representing the real and imaginay parts of a single bin. So, bin 0 (array indexes 0, 1) would represent the real and imaginary parts of the DC frequency. As you can see, bins 1 and 3 both have the same values, (except for the sign of the Im part), so I guess bin 3 is the first negative frequency, and finally indexes (4, 5) would be the last positive frequency bin.

Then to attenuate the frequency bin 1 this is what I do:

// 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;

For the actual tests I'm using a 1024-length FFT and I always provide all the samples so no 0-padding is needed.

// 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;
}

Obviously I'm doing something wrong but can't figure out what.

I don't want to use the FFT output as a means to generate a set of FIR coefficients since I'm trying to implement a very basic dynamic equalizer.

What's the correct way to filter in the frequency domain? what I'm missing?

Also, is it really needed to attenuate negative frequencies as well? I've seen an FFT implementation where neg. frequency values are zeroed before synthesis.

Thanks in advance.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

倦话 2024-09-10 20:05:18

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

这是图表：
IFFT of Product

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

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

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

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

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

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

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

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

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

There are two issues: the way you use the FFT, and the particular filter.

Filtering is traditionally implemented as convolution in the time domain. You're right that multiplying the spectra of the input and filter signals is equivalent. However, when you use the Discrete Fourier Transform (DFT) (implemented with a Fast Fourier Transform algorithm for speed), you actually calculate a sampled version of the true spectrum. This has lots of implications, but the one most relevant to filtering is the implication that the time domain signal is periodic.

Here's an example. Consider a sinusoidal input signal x with 1.5 cycles in the period, and a simple low pass filter h. In Matlab/Octave syntax:

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);

And here's the graph:
IFFT of product

The glitch at the beginning of the block is not what we expect at all. But if you consider fft(x), it makes sense. The Discrete Fourier Transform assumes the signal is periodic within the transform block. As far as the DFT knows, we asked for the transform of one period of this:
Aperiodic input to DFT

This leads to the first important consideration when filtering with DFTs: you are actually implementing circular convolution, not linear convolution. So the "glitch" in the first graph is not really a glitch when you consider the math. So then the question becomes: is there a way to work around the periodicity? The answer is yes: use overlap-save processing. Essentially, you calculate N-long products as above, but only keep N/2 points.

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

And here's the graph of ycorrect:

This picture makes sense - we expect a startup transient from the filter, then the result settles into the steady state sinusoidal response. Note that now x can be arbitrarily long. The limitation is Nproc > 2*min(length(x),length(h)).

Now onto the second issue: the particular filter. In your loop, you create a filter who's spectrum is essentially H = [0 (1:511)/512 1 (511:-1:1)/512]'; If you do hraw = real(ifft(H)); plot(hraw), you get:
hraw

It's hard to see, but there are a bunch of non-zero points at the far left edge of the graph, and then a bunch more at the far right edge. Using Octave's built-in freqz function to look at the frequency response we see (by doing freqz(hraw)):

The magnitude response has a lot of ripples from the high-pass envelope down to zero. Again, the periodicity inherent in the DFT is at work. As far as the DFT is concerned, hraw repeats over and over again. But if you take one period of hraw, as freqz does, its spectrum is quite different from the periodic version's.

So let's define a new signal: hrot = [hraw(513:end) ; hraw(1:512)]; We simply rotate the raw DFT output to make it continuous within the block. Now let's look at the frequency response using freqz(hrot):

Much better. The desired envelope is there, without all the ripples. Of course, the implementation isn't so simple now, you have to do a full complex multiply by fft(hrot) rather than just scaling each complex bin, but at least you'll get the right answer.

Note that for speed, you'd usually pre-calculate the DFT of the padded h, I left it alone in the loop to more easily compare with the original.

回复收藏 0 原文

浅浅淡淡 2024-09-10 20:05:18

您的主要问题是短时间间隔内的频率没有明确定义。对于低频尤其如此，这就是为什么您最容易注意到那里的问题。

因此，当您从声音序列中取出非常短的片段，然后对这些片段进行过滤时，过滤后的片段不会以产生连续波形的方式进行过滤，并且您会听到片段之间的跳跃，这就是您在此处生成喀哒声的原因。

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

首先我们将从 40Hz 的低通开始。因此，由于原始频率小于 40Hz，截止频率为 40Hz 的低通滤波器实际上不会产生任何影响，我们将得到几乎与输入完全匹配的输出。正确的？ 错误，错误，错误 - 这基本上是问题的核心。问题在于，对于较短的部分，27.5 Hz 的想法没有明确定义，并且无法在 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 没有太大关系。

回复收藏 0 原文

梦里兽 2024-09-10 20:05:18

您是否将直流频率样本的值衰减为零？在您的示例中，您似乎根本没有减弱它。由于您正在实现高通滤波器，因此您还需要将 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 返回时域，对每个输入音频缓冲区进行计算。这比您为每个样本所做的所有分支要高效得多。

Are you attenuating the value of the DC frequency sample to zero? It appears that you are not attenuating it at all in your example. Since you are implementing a high pass filter, you need to set the DC value to zero as well.

This would explain low frequency distortion. You would have a lot of ripple in the frequency response at low frequencies if that DC value is non-zero because of the large transition.

Here is an example in MATLAB/Octave to demonstrate what might be happening:

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');

Notice that in my code, I am creating an example of the DC value being non-zero, followed by an abrupt change to zero, and then a ramp up. I then take the IFFT to transform into the time domain. Then I perform a zero-padded fft (which is done automatically by MATLAB when you pass in an fft size bigger than the input signal) on that time-domain signal. The zero-padding in the time-domain results in interpolation in the frequency domain. Using this, we can see how the filter will respond between filter samples.

One of the most important things to remember is that even though you are setting filter response values at given frequencies by attenuating the outputs of the DFT, this guarantees nothing for frequencies occurring between sample points. This means the more abrupt your changes, the more overshoot and oscillation between samples will occur.

Now to answer your question on how this filtering should be done. There are a number of ways, but one of the easiest to implement and understand is the window design method. The problem with your current design is that the transition width is huge. Most of the time, you will want as quick of transitions as possible, with as little ripple as possible.

In the next code, I will create an ideal filter and display the response:

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');

Notice that there is a lot of oscillation caused by the abrupt changes.

The FFT or Discrete Fourier Transform is a sampled version of the Fourier Transform. The Fourier Transform is applied to a signal over the continuous range -infinity to infinity while the DFT is applied over a finite number of samples. This in effect results in a square windowing (truncation) in the time domain when using the DFT since we are only dealing with a finite number of samples. Unfortunately, the DFT of a square wave is a sinc type function (sin(x)/x).

The problem with having sharp transitions in your filter (quick jump from 0 to 1 in one sample) is that this has a very long response in the time domain, which is being truncated by a square window. So to help minimize this problem, we can multiply the time-domain signal by a more gradual window. If we multiply a hanning window by adding the line:

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

after taking the IFFT, we get this response:

So I would recommend trying to implement the window design method since it is fairly simple (there are better ways, but they get more complicated). Since you are implementing an equalizer, I assume you want to be able to change the attenuations on the fly, so I would suggest calculating and storing the filter in the frequency domain whenever there is a change in parameters, and then you can just apply it to each input audio buffer by taking the fft of the input buffer, multiplying by your frequency domain filter samples, and then performing the ifft to get back to the time domain. This will be a lot more efficient than all of the branching you are doing for each sample.

回复收藏 0 原文

辞取 2024-09-10 20:05:18

首先，关于标准化：这是一个已知（非）问题。 DFT/IDFT 需要每个因子中都有一个因子 1/sqrt(N)（除了标准余弦/正弦因子）（直接求逆）以使它们是对称的且真正可逆的。另一种可能性是将其中一个（正数或逆数）除以N，这是一个方便和品味的问题。 FFT 例程通常不执行此归一化，用户需要了解这一点并根据自己的喜好进行归一化。参见

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