DSP-通过FFT在频域中进行滤波 [英] DSP - Filtering in the frequency domain via FFT

问题描述

我一直在使用FFT的Exocortex实现进行一些尝试，但是我遇到了一些问题.

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

每当我在调用iFFT之前修改频率仓的幅度时，结果信号都会包含一些喀嗒声和pop声，尤其是当信号中存在低频(例如鼓或贝司)时.但是，如果我将所有垃圾箱衰减相同的因子，则不会发生这种情况.

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.

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

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

输出由成对的浮点数组成，每个浮点数代表单个容器的实部和虚部.因此，bin 0(数组索引0、1)将代表DC频率的实部和虚部.如您所见，bin 1和3都具有相同的值(Im部分的符号除外)，所以我想bin 3是第一个负频率，最后索引(4，5)将是最后一个正频率频率槽.

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.

然后衰减频率槽1，这就是我要做的:

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;

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

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.

我不想使用FFT输出来生成一组FIR系数，因为我试图实现一个非常基本的动态均衡器.

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?

还需要衰减负频率吗?我看过了FFT的实现，但否.频率值在合成之前归零.

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.

谢谢.

推荐答案

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

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

传统上，过滤是在时域中进行卷积实现的.没错，将输入信号和滤波器信号的频谱相乘是等效的.但是，当您使用离散傅立叶变换(DFT)(通过快速傅立叶变换算法实现速度)时，实际上是在计算真实频谱的采样版本.这有很多含义，但是与滤波最相关的一个含义是时域信号是周期性的.

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.

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

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

这是图形:

块开头的故障根本不是我们所期望的.但是，如果您考虑fft(x)，这是有道理的.离散傅立叶变换假定信号在变换块内是周期性的.据DFT所知，我们要求对此进行一个时期的转换:

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:

这导致在使用DFT进行过滤时的第一个重要考虑:您实际上是在实现圆形卷积 ，而不是线性卷积.因此，当您考虑数学时，第一个图形中的故障"并不是真正的故障.因此，问题就变成了:是否有一种方法可以解决周期性问题?答案是肯定的:使用重叠保存处理.本质上，您按上述方法计算N个长乘积，但仅保留N/2点.

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

这是ycorrect的图形:

这张图片很有意义-我们期望滤波器产生一个启动瞬态，然后结果稳定在稳态正弦响应中.请注意，现在x可以任意长.限制为Nproc > 2*min(length(x),length(h)).

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

现在讨论第二个问题:特定的过滤器.在循环中，创建一个滤波器，其频谱本质上是H = [0 (1:511)/512 1 (511:-1:1)/512]';如果执行hraw = real(ifft(H)); plot(hraw)，则会得到:

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:

很难看到，但是在图形的最左边缘有一堆非零点，然后在最右边缘有一堆非零点.使用Octave的内置freqz函数查看我们看到的频率响应(通过执行freqz(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)):

从高通包络线到零，幅度响应都有很多波动.同样，DFT固有的周期性正在起作用.就DFT而言，hraw一遍又一遍地重复.但是，如果您像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.

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

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

好多了.所需的包络在那里，没有所有的波纹.当然，现在的实现并不是那么简单，您必须通过fft(hrot)进行完全复杂的乘法运算，而不仅仅是缩放每个复杂的bin，但是至少您会得到正确的答案.

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.

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

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.

这篇关于DSP-通过FFT在频域中进行滤波的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！