用于数字滤波器设计的经典滤波器映射

信息处理 过滤器设计 有限脉冲响应 无限脉冲响应 巴特沃思 切比雪夫滤波器
2022-01-06 06:34:23

在四种经典模拟滤波器类型中:Butterworth、Chebyshev 、 EllipticBessel——与最小二乘法firls( ( ) 等?firpmremezmaxflat

(注意:这个问题不是关于简单 IIR 结构在环路滤波器、泄漏累加器、陷波滤波器中的常见和有用应用,也不是关于使用优化的 IIR 结构作为直接数字设计。我的问题是针对通过复制模拟经典来专门设计更高性能的低通、高通或带通结构——在我目前的狭隘观点之外,这样做可能会有实际用途)。

我被教导(弗雷德哈里斯和其他人)避免“复制模拟”的陷阱,因为那些经典类型的技术仅限于我们可以用相对较少数量的电感器和电容器来做的事情,而在数字世界中我们通过简单的延迟和乘法(以及用于多速率设计的非线性换向器,从而产生非常有效的 FIR 结构),拥有底层数学和可扩展性的全部功能。我对从 s 到 z 的映射的使用(在 1960 年代后期之前很常见,因为在模拟滤波器设计方面拥有丰富的知识)主要限于现有模拟滤波器的仿真和建模,而不是为创建新的数字滤波器常见的低通、高通和带通结构。

也就是说,除了建模和模拟之外,我可能会错过简洁且良好的实际应用,在这些应用中“复制模拟”会产生更好的解决方案。

最佳答案将列出常见低通、高通和带通滤波器设计的应用(不是 IIR 肯定会统治的陷波滤波器),特定于 FIR 滤波器的优化算法(包括优化的多速率结构)不可能超越性能,对于任何类过滤器类型(或证明为什么在这种情况下总是首选优化算法)。如果只是假设的陈述或怀疑,那么第二个最佳答案是至少展示这种映射滤波器的特定情况,以允许针对“优化的”直接数字解决方案进行测试。

3个回答

我相信,根据我们试图解决的问题,我们可以而且应该使用这两种方法:经典模拟滤波器设计的转换,以及使用优化方法在数字域中的直接设计。请注意,经典设计的属性非常有限:分段常数所需幅度响应、最小相位和定义为最大平坦度 (Butterworth) 或最小幅度近似误差 (Cauer = 椭圆) 或两者组合的最优性(两种类型的切比雪夫)。如果还需要其他任何东西,那么我们需要求助于其他方法。

但是,如果我们满足于仅具有(恒定增益)通带和阻带的简单设计,如果我们不关心相位响应,并且如果上述两个最优标准符合我们的目的,那么没有什么可以超越最低限度-源自经典模拟原型的相位 IIR 滤波器。为什么会这样?好吧,因为这些过滤器根据所选标准进行优化。根据定义,没有低阶过滤器也可以完成这项工作,因此,这样的过滤器可以实现最有效的实现。我们也不应该忘记设计非常简单,只是基于公式。IIR 滤波器设计的优化方法仍然存在问题,因为设计问题是高度非线性的,不能保证我们能找到全局最优解。我们总能找到一些解决方案——局部最优解——但这可能还不够好。

如果 - 正如您的问题所建议的那样 - 我们选择 FIR 滤波器作为所有滤波问题的解决方案,那么设计将比 IIR 滤波器的情况更容易,因为我们基本上只是在处理经过充分研究的多项式近似。但是,请注意,最有效的 FIR 滤波器设计方法会产生线性相位滤波器。如果需要高阶,这样的过滤器会引入巨大的延迟。在某些应用程序中,这种延迟可能令人望而却步。此外,如果需要高质量的幅度近似值(例如,几乎恒定的通带、非常高的阻带衰减),那么由此产生的高滤波器阶数也会导致内存和计算方面的大实现复杂性。

有一些用于非线性相位 FIR 滤波器设计的方法,但它们的应用范围较小,而且通常效率较低。此外,如果我们不关心相位,我们不想指定任何非线性所需的相位响应——这对于那些例程是必需的——但我们只想指定幅度响应并让相位自行处理. 在这种情况下,将需要最小相位设计。然而,即使有用于最小相位 FIR 滤波器设计的方法,这些方法效率也不高,并且在数值上不稳健,因此它们对于高滤波器阶数是不切实际的。

总而言之,在应用程序需要比经典模拟滤波器设计提供更多灵活性的所有情况下,无论如何我们都需要使用优化方法,而且我们通常会选择 FIR 滤波器,因为它们比(数字设计的)IIR 具有许多优势过滤器(例如,相对容易设计和实施,没有限制环,固有稳定性)。但如果问题可以通过经典设计解决,我认为没有理由使用其他任何东西,因为经典设计会产生可靠且非常有效的过滤器。

最后一点,我还没有提到贝塞尔过滤器。这些是我认为作为数字滤波器设计原型已过时的唯一。原因是由于双线性变换的频率扭曲,它们的最大平坦群延迟特性不会延续到数字域。如果需要最大平坦群延迟,我们需要直接在数字域中制定问题。Thiran 已经完成并解决了这个问题:

蒂兰,JP (1971-11-01)。“具有最大平坦群延迟的递归数字滤波器”。IEEE 电路理论汇刊。18 (6): 659–664。


例子:

我设计了一个15三阶椭圆低通滤波器,通带纹波为0.5dB 和 100 dB 的阻带衰减。椭圆滤波器的设计标准是与所需的恒定幅度响应的偏差最小。这对应于 Parks-McClellan 算法用于线性相位 FIR 滤波器的设计标准。事实证明,为了满足对幅度响应的相同要求,需要使用阶数线性相位滤波器667是必须的。两个滤波器的幅度响应如下所示:

在此处输入图像描述

最大通带和阻带近似误差几乎相同。但是请注意,出于实际目的,两个过滤器之间仍然存在一些明显的差异。首先,由于 IIR 滤波器的阻带纹波很少,因此在几个频带中的阻带衰减比 FIR 滤波器要大得多。其次,对于 IIR 滤波器,从通带到阻带的过渡更加清晰。这是极点相对接近单位圆的简单结果。当然,这样的极点在定点实现中会造成严重的问题。两个滤波器性能之间的另一个区别是线性相位 FIR 滤波器的预振铃,这在音频应用中是不可取的。

最后,线性相位 FIR 滤波器显然比最小相位 IIR 滤波器具有更大的延迟。它的群延迟是恒定的(667/2=333.5样本),而 IIR 滤波器的群延迟与频率有关。IIR 滤波器的平均通带群延迟约为20样品(见下图)。如果相位的线性不是必需的,我们不妨使用与线性相位 FIR 滤波器具有相同幅度的最小相位 FIR 滤波器。我在这个答案中解释了次优最小相位滤波器的计算如下图所示,最小相位 FIR 滤波器的群延迟与 IIR 滤波器的群延迟非常相似。

关于复杂性,很明显15三阶 IIR 滤波器比相应的阶 FIR 滤波器实现起来要便宜得多667.

在此处输入图像描述

和模拟,其中“复制模拟”会产生更好的解决方案。

这有点漏掉了重点。并不是人们很关心匹配或复制“模拟”,而是数字 IIR 滤波器具有一些非常好的和有用的特性。

例如在音频中,IIR 滤波器非常常见,我每天都使用 Butterworths。

最直接的原因只是过滤器长度。当采样率为 48 kHz 时,要在 40 Hz 下做任何有意义的事情,FIR 滤波器需要有数千个抽头长度。在某些情况下,您可以使用多速率过滤器来非常有效地实现这一点,但我不确定这是否普遍适用。另一种可能的选择是基于 FFT 的过滤器(重叠相加等),但这在延迟和效率之间提出了一个重要的权衡。

如果你坚持使用线性相位滤波器,延迟很快就会变得非常大,并且你会失去因果关系。从技术上讲,您可以通过简单地反转单位圆外的零点来将任何线性相位 FIR 转换为最小相位 FIR,但这在高阶数字上是一个尴尬的过程。

人耳对单耳相位失真(例如最小相位)也相当不敏感,只要其原因。然而,它对预振铃非常敏感,因此线性相位滤波器在感知上可能存在问题。

一个更“哲学的原因”是人类听觉或多或少地使用对数频率轴,而 IIR 滤波器更适合于此。例如,8kHz 的 IIR 倍频程带通滤波器看起来与 125Hz 的滤波器完全相同(中心频率除外)。FIR 滤波器将大不相同。

IIR 非常适合分频:奇阶巴特沃斯低通和高通滤波器总和为平坦的频率响应(尽管它通常是全通的)。对于偶数订单,您可以使用 Linkwitz Riley 过滤器(基于 Butterworths)。您可以改用线性相位 FIR,但同样会遇到延迟和因果关系问题。

像实时可调高通和低通滤波器这样的操作也很容易:动态设计巴特沃斯滤波器或简单地列出极点位置并进行插值非常容易。过滤器顺序始终保持不变,并且保证您具有相同的斜率和形状。

换句话说:打开你家中的任何扬声器,你会听到很多巴特沃思和朋友在行动:-)

示例 1

这只是一个非常简单的高通,您可以在任何花园品种有源扬声器中使用。以 48 kHz 采样的 40 Hz。我使用firls(). 匹配是手动完成的,所以它们看起来一样。我需要在 FIR 滤波器上进行大约 6000 个抽头才能使其达到预期效果

% IIR
[z,p,k] = butter(6,40*2/48000,'high');
%FIR
h=firls(6000,[0 10 65 fs/2]*2/fs,[0 0 1 1])';

传递函数看起来像这样。 在此处输入图像描述

对我来说,FIR 在延迟、内存占用方面似乎非常逊色。一个 6 阶 BW 只有 13 个系数和 6-8 个状态变量。您可以通过利用零都在z=1.

示例 2:

这是对前一个示例的扩展,它使高通的截止频率可以实时调整。这会创建一个“滑动高通”滤波器,这在密封的智能扬声器中很常见。我试图让需求变得现实,就像你在真实产品中发现的那样。

Butterworth 允许非常高效的实现,没有可听伪影、没有固有延迟、非常低的内存占用和非常低的 CPU 消耗。

%% Script to implement a sliding high pass filter, that can be adjusted on the fly
% This type of "sliding highpass" is typically used in smart speakers with a
% "closed box" topology to control trasnducer excursion at high output
% volume
%
% Requirements:
%  - sample rate: 48 kHz
%  - frame size: 128 samples
%  - cutoff frequency varies from 40 Hz to 100 Hz
%  - highpass slope: 36 dB/octave. In other words level at half the cutoff
%      frequency should be <= -36 dB
%  - Level at cutoff frequency not less than -3dB
%  - passband ripple: < .01 dB above 200 Hz (where spectral perception is
%                      more sensitive)
%  - cutoff frequency is updated once per frame
%  - no additional latency
%  - smooth updating filter, no pops or clicks  
%
% Implementation
%  - 6th order butterworth
%  - pole location and filter gain as a function of frequency are 
%    approximated as a polynomial fit. At these low frequencies, a linear
%    fit (1rst order polynomial) works perfectly.
%  - The filter coefficients are calcuated once per frame based on the
%    current cutoff frequency
%  - Filter is implemented as cascaded second order sections in Direct Form
%    I. For Direct Form I, the state variable are simply the inputs and
%    outputs, so updateing the filter deosn't create a discontinuity in the
%    state variables.
%
%  Test
%  - input signal is a 50 Hz sine wave, low frequency sine waves are very
%    sensitive to artifacts
%  - cutoff frequency input is a mixture of step function, up and down sweeps
%    and uniformly distributed random numbers between 40Hz and 100Hz
%


%% Table up the poles and the gain of the filters, perform polynomial fit
ord = 6;
fs = 48000;
fr = (40:.1:100)';
nfr = length(fr);
p0 = cell(nfr,1);  % pole locations
k0 = zeros(nfr,1); % filter gains
for i = 1:nfr
  [z,p,k] = butter(ord,2*fr(i)/fs,'high');
  p0{i} = p(imag(p)>0);
  k0(i) = k;
end
p0 = [p0{:}].'; % convert cell array to regular array

%% do a simple linear for poles and gains
% we fit the real part and the gain in (1-x) since they are very close to 1

ppPoles = cell(3,2);
for i = 1:3
  ppPoles{i,1} = polyfit(fr,1-real(p0(:,i)),1); % 1 minus real part
  ppPoles{i,2} = polyfit(fr,imag(p0(:,i)),1); % imginary part
end

% and the gains in 1-k
ppGain = polyfit(fr,1-k0,2);

%% test the polyfit, calculate the resulting transfer functions at a few
% frequencies
frTest = 40:10:100;
nfrTest = length(frTest);
nFFT = 16384*2;
d0 = zeros(nFFT,1); d0(1) = 1;
fy = zeros(nFFT,nfrTest);
sos = zeros(ord/2,6);
sos(:,1) = 1; sos(:,2) = -2; sos(:,3) = 1;
for i = 1:nfrTest
  f =  frTest(i);
  for ip = 1:3
    x = 1-polyval(ppPoles{ip,1},f); % real part
    y = polyval(ppPoles{ip,2},f); % imginary part
    sos(ip,4:6) = [1 -2*x x^2+y^2];
  end
  k = 1-polyval(ppGain,f);

  fy(:,i) = fft(k*sosfilt(sos,d0));
end
% this checks out fine

%% Real time part, preperation
frameSize = 128;  % frame size
frameRate = fs/frameSize; % number of frames per second

% let's do 10 seconds but an integer number of frames
nx = 10*fs;
numFrames = floor(nx/frameSize);
nx = numFrames*frameSize;

% build a test signal: 50 hz sine wave
xin = sqrt(.5)*sin(2*pi*(0:nx-1)'*50/fs); % 50 Hz sine wave at -3 dB
% build control frequency signal, one frequency per frame
frInput = 40*ones(numFrames,1);
% let's do a bit of rectangular switching plus some random stuff
t = (1:frameRate);
frInput(frameRate+t) = 100;
frInput(3*frameRate+t) = 70;
% up and down sweep
frInput(5*frameRate+t) = linspace(40,100,frameRate);
frInput(6*frameRate+t) = flip(linspace(40,100,frameRate));
% some random numbers for good measure
frInput(7*frameRate+1:end) = 40+60*rand(3*frameRate,1);
% 
% in order to smooth the very aprupt frequency transitions in the test
% vector, we smooth the frequency input with a time constant to 30ms
timeConstant = 0.03;
frSmooth = exp(-1./(timeConstant*frameRate));
frCurrent =40;

%% Real time over all frames
% we implement this as driect form I so the states are always guaranteed to
% be continous
signalState = zeros(2,4);  % filter state, we need total of 8
freqState = frInput(1); % cutoff frequency states
xout = 0*xin; % initialize output
t = 1:frameSize; % time vector
t0 = 0;  % current time
yy = [xout xout xout];
for iFrame = 1:numFrames
  % get frequency input and apply smoothing
  frCurrent = frSmooth*frCurrent + (1-frSmooth)*frInput(iFrame);
  % calculate filter gain, grab input and scale it
  k = 1-polyval(ppGain,frCurrent);
  y = k*xin(t0+t);

  % over all biquads
  x1 = signalState(1,1); % grab input state for the first biquad
  x2 = signalState(2,1);
  for iPole = 1:3
    % calculate the filter coeefficents
    pReal = 1-polyval(ppPoles{iPole,1},frCurrent); % real part
    pImag = polyval(ppPoles{iPole,2},frCurrent); % imaginary part
    a1 = -2*pReal;                % filter coefficient "a1"
    a2 = pReal*pReal+pImag*pImag; % biquad coefficient "a2"
    % grab the output state
    y1 = signalState(1,iPole+1);
    y2 = signalState(2,iPole+1);
    % inner loop. Here efficieny is the most important
    for i = t
      x0 = y(i);  % get input sample
      y0 = x0-2*x1+x2-a1*y1-a2*y2; % DF1 Butterworth stage
      % update state
      x2 = x1; x1 = x0;
      y2 = y1; y1 = y0;
      y(i) = y0; % write output
    end % end sample loop
    % save the output state of this stage
    signalState(1,iPole) = x1;
    signalState(2,iPole) = x2;
    % grab the input state for the next stage
    x1 = signalState(1,iPole+1);
    x2 = signalState(2,iPole+1);
    % now write the output state
    signalState(1,iPole+1) = y1;
    signalState(2,iPole+1) = y2;

    yy(t0+t,iPole) = y;
      
  end % end poles/stages
  % write output
  xout(t0+t) = y;
  t0 = t0 + frameSize;
  
end % end frame

plot(xout);

如何在 FIR 滤波器方法中大量减少资源需求

Matt 和 Hilmar 提供的两个答案都非常出色,并且在回答问题时提供了很好的洞察力。我赞成 Hilmar 的回答是正确的,因为它既能很好地触及并很好地展示了要点,但我对他提供的效率声明尚不信服。在选择答案之前,需要解决这个问题(我的错误或删除的那些声明)。

Hilmar 滤波器的 FIR,每个时钟周期有 8.45 个乘法器

我在下面显示的内容由于他选择的信号带宽和紧密的过渡宽度相结合,我无法使用 Matt 的滤波器(如果没有可以完成的等效技巧,我不会感到惊讶,但我还没有看到它们)。但是,我能够显着减少 Hilmar 作为最小二乘 FIR 滤波器的实现。Hilmar 强调了许多其他优点,使他的实现对某些应用程序具有吸引力,但我确实想澄清 FIR 滤波器在每个时钟周期需要数千个乘法器的误解(或者我自己的误解/错误)。原因是多速率信号处理——以及通过以尽可能低的速率运行所需的信号处理的选定部分可以获得的优势。

过滤性能比较

首先,我将重点介绍 Hilmar 的 IIR 滤波器与我的 FIR 滤波器实现的比较结果。Hilmar 没有详细说明要求,但我使用了 100 dB 的 16 位 48 KHz 音频动态范围(阻带与 Matt 在他的帖子中使用的一致)。

下面显示了过渡带的对数频率图,与 Hilmar 帖子中的图一致,除了我的 IIR 和 FIR 实现版本。

40 Hz 高通比较

下面是一个类似于 Matt 的演示文稿的四边形图表,但同样适用于我的 Hilmar 的 IIR 和 FIR 滤波器版本。对于这两种滤波器实现,我使用 28 位系数量化来包含可能的量化效果,并发现两者都需要 28 位,以便在相同浮点结果显示的这个比例上不偏离:

四元图过滤器比较

6 阶 40 Hz HPF IIR 实现

从我认为相当于 Hilmar 的实现开始;下面是我为 6 阶 Butterworth IIR 滤波器确定的 2 阶双二阶系数,量化为 28 位并归一化以将所需的乘法次数减少到 9(总共 11 个,但两个只需要位移)。值得注意的是,所有这些乘法都以 48 KHz 运行!

[[ 0.98993582 -1.97987163  0.98993582  1.         -1.98990852  0.98993579]
 [ 1.         -2.          1.          1.         -1.99259523  0.99262254]
 [ 1.         -2.          1.          1.         -1.99726595  0.99729332]]

下面是一个框图,显示了这种 IIR 方法的实现,以防可以确定进一步的简化。我对 FPGA 实现的偏好是使用所示的这种转置的直接形式 2 表示,因为它避免了加法器树,并允许在遇到时序延迟问题之前以可能的最高速率运行实现(但这在 48KHz 音频滤波器中肯定不是问题)。

六阶 IIR

5533 阶 40 Hz HPF FIR 实现

由于抽取和插值技术的强大功能,我们可以以尽可能低的速率运行内核滤波操作,从而显着减少每个时钟周期的乘法器操作 - 从而最大限度地减少资源,并显着降低功耗(与时钟频率;P=CV2f)。当然,我们放弃的是我将详细说明的内存需求,以及总体延迟,出于这些原因,我们可能会支持 Hilmar 建议的 IIR 方法,但看起来我们实际上并不需要在每个时钟周期上使用 1000 次乘法器。下面的框图的要点是,通过创造性的抽取和插值技术,我们将大量的处理转移到了 1200 Hz 而不是 48KHz!每个 48 KHz 时钟周期所需的实数 28 x 16 乘法的总数以红色计算。

带抽取和插值的 40 Hz FIR

上面的框图在功能上逐位等效于以 48 KHz 运行的 5534 抽头线性相位 FIR 滤波器,每个时钟周期需要 2217 次乘法:

FIR 5534 水龙头

我在这里所做的是实现最小二乘线性相位低通这是滤波器的群延迟(57 ms)。因此,如果这种延迟是一个问题,或者是保存 2.7K 样本所需的内存,那么 Hilmar 提出的 IIR 将是一个令人信服的选择。

脉冲响应

下面显示了两种实现的脉冲响应,在幅度的对数尺度上突出了差异。在这里,我们看到了线性相位 FIR 滤波器的可怕“主要回声”,我从 RBJ、Hilmar 和其他发烧友那里了解到,这在音频世界中是纯粹的邪恶。请注意,在 Matt 和 Hilmar 的实现中,IIR 滤波器方法都没有前导回波并且总体延迟显着减少。

脉冲响应

这是抽取/插值方法的第一次切割,以显示 FIR 与 IIR 滤波器的真实比较。通过不同的组合或进一步降低时钟频率,资源需求可能会进一步降低(当然,我们可以在 200 甚至更低的频率下运行内部 FIR 以获得 40 Hz 的高通!)。我们也可以使用带有抽取/插值的最小相位 FIR 来做到这一点,但是我们失去了能够通过简单地减去延迟将低通转换为高通的好特性(在 Hilmar 的情况下)。匹配减法延迟的工作可能非常具有挑战性,而且如果我们要这样做,那么我们可以使用 IIR LPF 来做到这一点,并提供相同的抽取/插值优势。这是一个线性相位 FIR,这一点对于这个例子很重要。

资源需求摘要

FIR 的资源需求摘要如下:

  • 每 48 KHz 时钟周期 8.45 次实数乘法

  • 每 48 KHz 时钟周期 13.35 个实数加法

  • 2747 x 16 位 RAM 用于输入信号

  • (3x8 + 3x5 + (135-1)/2+1 ) = 104 x 24 位 RAM 用于系数

  • 207 x 16 位 RAM 用于过滤器状态存储

因此,对于内存和 57 ms 延迟不是问题,并且功耗和乘法器资源最受关注的情况 - FIR 方法提供了一个非常引人注目的选项,即使对于 Hilmer 建议的这种具有挑战性的实现(动态可调性的功能也很容易完成)假设它完全被使用的核心 135 抽头 FIR 滤波器的响应所覆盖)。这是初稿,如果需要,我相信资源需求可以进一步显着降低,因为我们仍在以相对较高的 1.2 KHz 速率运行具有 40 Hz 截止频率的滤波器。可能仍然是 FIR 方法无法轻易消除但最小化的问题是主要回波。这一切都非常有趣,我非常感谢 Matt 和 Hilmar 的贡献和想法。

有关诸如此类的多速率技术的完整详细信息,请参阅 fred harris 最近出版的多速率信号处理第二版(我现在找不到任何低于 500 美元的第一版副本)。