IIR Bandpass Filters Using Cascaded Biquads

翻译 ALLEN ⋅ 于 2019-08-02 22:31:25 ⋅ 最后回复由 ALLEN 2020-06-29 10:44:13 ⋅ 68 阅读

在之前的文章 [1] (参考文献 [1])中,我们实现了使用二阶 IIR 过滤器级联的低通 IIR 滤波器,即 biquad。本文提供了一个 Matlab 函数,可以对 Butterworth 带通 IIR 滤波器执行相同的操作。与传统方法相比,基于 biquad 的带通滤波器对系数量化 [2](参考文献 [2]) 的灵敏度较低。这在设计窄带滤波器时变得很重要。

图1 显示了使用直接形式 II 结构[3,4]的 biquad 切片框图。它是由两个反馈系数 a,三个前馈系数 b 和增益 k 定义的。函数 biquad_bp (在附录中列出)根据低通原型、中心频率、带宽和采样频率的顺序计算 biquad系数和增益。对于一个 N 阶的低通原型,带通滤波器的阶数为 2N,需要N 个四分体。
file
图 1 Biquad 切片框图

例 1.

下面是一个基于三阶低通原型的带通滤波器的函数调用示例和结果:

    N= 3;               % order of prototype LPF (number of biquads in BPF)
    fcenter= 20;        % Hz center frequency
    bw= 4;              % Hz -3 dB bandwidth
    fs= 100;            % Hz sample frequency
    [a,K]= biquad_bp(N,fcenter,bw,fs)
        a = 1.0000  -0.3835  0.8786
            1.0000  -0.5531  0.7757
            1.0000  -0.7761  0.8865
                                                                 .        
        K = 0.1254  0.1122  0.1114

该函数在矩阵 a 中生成 N = 3 组 biquad 反馈系数,并为每个 biquad 段增加一个增益 K。计算增益 K 使每个部分的响应幅度为1.0,大约为$f_{center}$。我们可以用 a 和 K 求出四轴飞行器的频率响应。我们稍后将展示,biquad 的前馈系数都是相同的:b = [1 0 -1]。我们将每个biquad 的参数赋值如下:

    b= [1 0 -1];            % feed-forward coeffs of each biquad
    a_bq1= a(1,:);          % biquad 1 feedback coeffs = [1 -.3835 .8786]
    a_bq2= a(2,:);          % biquad 2 feedback coeffs = [1 -.5531 .7757]
    a_bq3= a(3,:);          % biquad 3 feedback coeffs = [1 -.7761 .8865]
    K1= K(1); K2= K(2); K3= K(3);    % biquad gains

现在我们可以计算每个 biquad 的频率响应。

    [h1,f]= freqz(K1*b,a_bq1,2048,fs);    % frequency response
    [h2,f]= freqz(K2*b,a_bq2,2048,fs);
    [h3,f]= freqz(K3*b,a_bq3,2048,fs);
    H1= 20*log10(abs(h1));                % dB-magnitude response
    H2= 20*log10(abs(h2));
    H3= 20*log10(abs(h3));

这些响应如图2所示。计算级联 biquads 的响应:

    h= h1.*h2.*h3;
    H= 20*log10(abs(h));

整个响应如图3所示。给定前馈系数 b = [1 0 -1],各 biquad 的框图如图4所示。系数 a0 = 1.0 没有使用。注意,反馈增益有负号,所以给定 a1 = -0.3835 和 a2 = 0.8786,反馈增益分别为 0.3835 和 -0.8786。再次查看图2,注意第一个和最后一个 biquad 在某些频率上的响应幅度大于 1.0 (0 dB),即使在$f_{center} = 20$ Hz附近的响应是 0 dB。在实现定点滤波器时,我们需要考虑这一点以避免溢出。在本例中,如果输入在第一个 biquad 的响应峰值(约在 22 Hz)有一个较大的分量,则最有可能发生溢出。为了防止溢出,我们可以在第一个 biquad 允许额外的空间;将$2^{nd}$ biquad与第一个 biquad 交换;或降低增益$K_{1}$,同时保持 biquad 增益的乘积不变。
file
图2 在 N = 3,中心频率 = 20 Hz,带宽 = 4 Hz,fs = 100 Hz时,各biquad 段的响应
file
图3 带通滤波器的整体响应
file
图4 带通滤波器的 Biquad 部分

BIQUAD_BP 是怎样工作的

在之前提交的IIR 带通滤波器 [5] 文章 中,我讲述了零极点的带通响应形式为:

$$H(z)=K\frac{(z+1)^N(z-1)^N}{(z-p_1)(z-p_2)...(z-p_{2N})}\qquad (1)$$

其中 N 是 原型低通滤波器的阶数,文章中描述了 $p_{k}$ 的计算过程。对于一个带通滤波器,2N 个极点作为复数共轭对的形式出现。比如,例 1 的滤波器的极点如图 5 所示。在 z = -1 和 z = +1 也各自存在着 N 个零点。我们的目标是把 H(z) 转换为多个二阶部分(biquads)的级联。通过把 H(z) 写成 具有复数共轭极点对的N 个部分的成绩来实现这个目标:

$$H(z)=K_1\frac{(z+1)(z-1)}{(z-p_1)(z-p_1^*)}\cdot K_2\frac{(z+1)(z-1)}{(z-p_2)(z-p_2^*)}\cdot…\cdot K_N\frac{(z+1)(z-1)}{(z-p_N)(z-p_N^*)}\qquad (2)$$

其中 $p^{*}_{k}$ 是极点 $p_{k}$ 的复数共轭。注意到,我们给每个 biquad 在 z = -1 和 z = +1 各自赋了一个零。我们可以在不同的在值上分配零点,但是上述方式是比较直接的,并且给出了合理的结果。对第 k 个 biquad 部分的分子分母进行展开,我们得到:

$$H_k(z)=K_k\frac{z^2-1}{z^2-(p_k+p_k^*)z+p_kp_k^*} =K_k\frac{z^2-1}{z^2+a_1z+a_2}$$

其中 $a_{1}= -2*real(p_{k})$ 和 $a_{2}= |p_{k}|^2$。分子和分母同除以$z^2$,我们得到:

$$H_k(z)=K_k\frac{1-z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}\qquad (3)$$

因为我们对每个 biquad 分配了相同的零点,分子(前馈)系数 b = [1 0 -1] 对所有 N 个biquad是相同的。 对这些 biquad 系数进行相加,我们有:

$$b = [1 0 -1]$$$$a = [1 -2*real(p_{k}) |p_{k}|^2] $$

为了得到增益$K_{k}$,我们让每个 biquad 在滤波器几何平均频率$f_{0}$处具有增益 1 。为实现这个,我们在$f_{0}$处计算$H_{k}(z)$,并且设置$K_{k} = 1/|H(f_{0})|$。为了得到$f_{0}$,我们定义$f_{1}= f_{center} – bw/2$和$f_{2}= f_{center}+bw/2$。然后,$f_{0}=\sqrt{f_{1}*f_{2}}$。注意,对于窄带滤波器,$f_{0}$位于$f_{center}$附近。

计算$K_{k}$的 Matlab 代码如下:

    f0= sqrt(f1*f2);
    h= freqz(b,a,[f0 f0],fs);     % frequency response at f = f0
    K= 1/abs(h(1));

具有系数 b,a 和 增益 K 的 biquad 部分如图 4 所示。以下是使用 biquad_bp (前4个步骤的详细信息在文献 [5])滤波器合成步骤的步骤:

  1. 获得 Ωc = 1 rad/s 的低通模拟原型滤波器的极点。
  2. 给定上下 -3 dB 频率的数字带通滤波器,获得模拟带通滤波器(pre-warping)的相应频率。
  3. 把模拟低通极点转化为模拟带通极点。
  4. 使用双线性变换,把 s 平面的极点转换为 z 平面的极点。
  5. 使用上面的方程 4 计算每个 biquad 的反馈(分母)系数。
  6. 计算每个 biquad 的增益 K。
    file
    图 5. 例1的带通滤波器的z平面极点和零点

    例2-带量化系数的窄带滤波器

    让我们试着实现一个带量化的分母系数的窄带通滤波器。我们会使用一个采样频率为100Hz的 0.5 带宽。 使用$3^{rd}$阶和中心频率为22Hz的低通原型,函数调用和结果为:

    N= 3;                % order of prototype LPF
    fcenter= 22;         % Hz center frequency
    bw= .5;              % Hz -3 dB bandwidth
    fs= 100;             % Hz sample frequency
    [a,K]= biquad_bp(N,fcenter,bw,fs)
     a = 1.0000  -0.3453  0.9844
         1.0000  -0.3690  0.9691
         1.0000  -0.3983  0.9845
                                                          .
     K = 0.0157  0.0155  0.0155

    现在给分子系数赋值,并且用12位数对分母系数进行量化:

     b= [1 0 -1];                        % numerator coeffs
     nbits= 12;                          % bits per unit amplitude
     a_quant= round(a*2^nbits)/2^nbits;  % quantize denominator coeffs

    例 1 计算了频率响应,其响应如图 6 所示。图7展示了 12 位和 8 位分母系数时的1 dB/division的响应。就如你所见的,我们使用 12 位系数得到好的结果以及使用8位系数得到的一些反应降解。图 8 展示了三个具有六个极点的过滤器,相对于例 1(图 5)的少窄带滤波器,这些极点非常接近于单位圆。由于极点落在单位圆外,带量化系数的传统带通实现极易变得不稳定。
    file
    图 6. N = 3, bw = 0.5 Hz,以及分母系数 = 12 位的带通响应。

file
图 7. 带通滤波器响应 vs. 分母系数量化.

file
图 8. 例 2 三个具有 6 个 z 平面极点的滤波器。

APPENDIX MATLAB FUNCTION BIQUAD_BP

本程序中的一些公式是在之前的带通 IIR 滤波器上开发出来的[5]。本程序按原样提供,不作任何保证或保证。对于因使用或误用程序而造成的任何损坏或损失,作者不承担任何责任。

%function [b,a]= biquad_bp(N,fcenter,bw,fs)  4/7/19  Neil Robertson
% Synthesize IIR Butterworth Bandpass Filters using biquads
%
% N= order of prototype LPF = number of biquads in BPF
% fcenter= center frequency, Hz
% bw= -3 dB bandwidth, Hz
% fs= sample frequency, Hz
% a= matrix of denominator coefficients of biquads.
%    each row contains the denom coeffs of a biquad.
%    There are N rows.
% K= vector of gain blocks for each biquad.Length = N.
%
% Note:numerator coeffs of each biquad = K(k)*[1 0 -1]
%
function [a,K]= biquad_bp(N,fcenter,bw,fs)
f1= fcenter- bw/2;                % Hz lower -3 dB frequency
f2= fcenter+ bw/2;                % Hz upper -3 dB frequency
if f2>=fs/2;
    error('fcenter+ bw/2 must be less than fs/2')
end
if f1<=0
    error('fcenter- bw/2 must be greater than 0')
end
% find poles of butterworth lpf with Wc = 1 rad/s
k= 1:N;
theta= (2*k -1)*pi/(2*N);
p_lp= -sin(theta) + j*cos(theta);
% pre-warp f0, f1, and f2 (uppercase == continuous frequency variables)
F1= fs/pi * tan(pi*f1/fs);
F2= fs/pi * tan(pi*f2/fs);
BW_hz= F2-F1;                % Hz -3 dB bandwidth
F0= sqrt(F1*F2);             % Hz geometric mean frequency
% transform poles for bpf centered at W0
% pa contains N poles of the total 2N -- the other N poles not computed
% are conjugates of these.
for i= 1:N;
    alpha= BW_hz/F0 * 1/2*p_lp(i);
    beta= sqrt(1- (BW_hz/F0*p_lp(i)/2).^2);
    pa(i)= 2*pi*F0*(alpha +j*beta);
end
% find poles of digital filter
p= (1 + pa/(2*fs))./(1 - pa/(2*fs));    % bilinear transform
% denominator coeffs
for k= 1:N;
    a1= -2*real(p(k));
    a2= abs(p(k))^2;
    a(k,:)= [1 a1 a2];                 % denominator coeffs of biquad k
end
b= [1 0 -1];                           % biquad numerator coeffs
% compute biquad gain K for response amplitude = 1.0 at f0
f0= sqrt(f1*f2)                        % geometric mean frequency
for k= 1:N;
    aa= a(k,:);
    h= freqz(b,aa,[f0 f0],fs);         % frequency response at f = f0
    K(k)= 1/abs(h(1));
end

参考文献

Robertson, Neil, “Design IIR Filters Using Cascaded Biquads”, Feb, 2018, https://www.dsprelated.com/showarticle/1137.php
Oppenheim, Alan V. and Shafer, Ronald W., Discrete-Time Signal Processing, Prentice Hall, 1989, section 6.9.
Sanjit K. Mitra, Digital Signal Processing, 2nd Ed., McGraw-Hill, 2001, section 6.4.1
“Digital Biquad Filter”, https://en.wikipedia.org/wiki/Digital_biquad_filter
Robertson, Neil, “Design IIR Bandpass Filters”, Jan, 2018,
https://www.dsprelated.com/showarticle/1128.php


原文地址:https://www.dsprelated.com/showarticle/1257.php...

译文地址:...

追求梦想,做最好的自己

本帖由 ALLEN 于 10个月前 置顶
回复数量: 2
  • zhangzs
    在线: 18350 分钟
    zhangzs
    2019-09-10 09:44:13

    这个文章是从哪弄得?

  • ALLEN
    在线: 113141 分钟
    ALLEN MOD 主攻数据挖掘,机器学习,数字信号处理,语音信号处理
    2019-09-10 09:46:20

    在dsprelated.com

暂无评论~~
  • 请注意单词拼写,以及中英文排版,参考此页
  • 支持 Markdown 格式, **粗体**、~~删除线~~、`单行代码`, 更多语法请见这里 Markdown 语法
  • 支持表情,使用方法请见 Emoji 自动补全来咯,可用的 Emoji 请见 :metal: :point_right: Emoji 列表 :star: :sparkles:
  • 上传图片, 支持拖拽和剪切板黏贴上传, 格式限制 - jpg, png, gif
  • 发布框支持本地存储功能,会在内容变更时保存,「提交」按钮点击时清空
  请勿发布不友善或者负能量的内容。与人为善,比聪明更重要!
Ctrl+Enter