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的实部和虚部.因此,bin 0(数组索引0,1)将代表DC频率的实部和虚部.正如你所看到的,第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系数的方法,因为我正在尝试实现一个非常基本的动态均衡器.

在频域中过滤的正确方法是什么？我错过了什么？

此外,是否真的需要衰减负频率？我已经看到了一个带有neg的FFT实现.频率值在合成之前归零.

提前致谢.

Answer 1

有两个问题:使用FFT的方式和特定的过滤器.

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

这是一个例子.考虑一个x周期为1.5个周期的正弦输入信号,以及一个简单的低通滤波器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进行过滤时,这导致了第一个重要的考虑因素:实际上是在实现循环卷积,而不是线性卷积.因此,当您考虑数学时,第一个图中的"故障"并不是真正的故障.那么问题就变成了:有没有办法解决周期问题？答案是肯定的:使用重叠保存处理.基本上,您如上所述计算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

这是以下图表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,但至少你会得到正确的答案.

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

Answer 2

您的主要问题是频率在短时间间隔内没有明确定义.对于低频率尤其如此,这就是为什么你注意到那里最常见的问题.

因此,当您从声音列中取出非常短的片段,然后对其进行过滤时,过滤后的片段将不会以产生连续波形的方式进行过滤,并且您会听到片段之间的跳跃,这就是您在此处产生的点击.

例如,取一些合理的数字:我开始时波形为27.5 Hz(钢琴上的A0),数字化为44100 Hz,它看起来像这样(红色部分为1024个样本长):

替代文字http://i48.tinypic.com/zim802.png

所以首先我们将从40Hz的低通开始.因此,由于原始频率低于40Hz,具有40Hz截止频率的低通滤波器实际上不会产生任何影响,我们将获得几乎与输入完全匹配的输出.对？ 错,错,错 - 这基本上是你问题的核心.问题在于,对于短截面,27.5Hz 的想法没有明确定义,并且在DFT中无法很好地表示.

通过查看下图中的DFT,可以看出27.5 Hz在短段中没有特别意义.请注意,虽然较长的段的DFT(黑点)在27.5 Hz处显示峰值,但较短的点(红点)不会.

alt text http://i50.tinypic.com/14w6luw.png

显然,然后在低于40Hz的频率下滤波将仅捕获DC偏移,并且40Hz低通滤波器的结果在下面以绿色显示.

替代文字http://i48.tinypic.com/2vao21w.png

蓝色曲线(以200 Hz截止频率拍摄)开始更好地匹配.但请注意,不是低频使其匹配良好,而是包含高频.直到我们在短段中包含可能的每个频率,高达22KHz,我们才能最终得到原始正弦波的良好表示.

所有这一切的原因是27.5 Hz正弦波的一小段不是 27.5 Hz的正弦波,它的DFT与27.5 Hz没有太大关系.