一文读懂线性相位滤波器和最小相位滤波器

前一篇博客 《一文读懂滤波器的线性相位，全通滤波器，群延迟》 详细解释了线性相位滤波器，与之对应的是最小相位滤波器，本文就最线性相位滤波器和最小相位滤波器进行比较。并提出问题： 线性相位滤波器一定更适合音频吗？ 先说结论：

1. 线性相位滤波器

• 优点：能实现线性相位
• 缺点：pre-ringing效应，听觉上类似chirp的效果。并且具有较高的延迟。

2. 最小相位滤波器

• 优点：延时低，相对于IIR对相位影响更小。
• 缺点：会带来相位的变化，只能用来补偿幅频响应，不能补偿相频响应。

1. 线性相位和最小相位滤波器

对于递归实现的滤波器（即IIR），不存在线性相位，因为递归滤波器无法生成对称的脉冲响应。然而对于离线信号，可以通过两次使用递归滤波器来实现零相位。例如，对于离线存在的信号$x(\cdot )$，可以分别向前和向后应用递归滤波器。这样做得到了滤波器幅度响应的平方，并将相位响应归零。通过以下公式可以看出，零相位滤波器（zp）和最小相位滤波器（mp）幅度响应相同，而零相位滤波器相位响应为零。零相位滤波器是线性相位滤波器的一种，其每个频点的群延迟都是零。
Z { f l i p ( x ) } ≜ X ( z − 1 ) Z\{flip(x)\} \triangleq X(z^{-1}) Z{flip(x)}≜X(z−1)
$$Y_{zp}(z)=H(z^{-1})[H(z)X(z)]$$
$$Y_{mp}(z)=H(z)[H(z)X(z)]$$

1.1 椭圆滤波器

选择椭圆滤波器是因为它在截止频率附近具有高度非线性的相位响应，从而导致额外的延迟，在该频率下可以听到“振铃”（ringing）。截止频率选择 2kHz，是因为人耳敏感的频率。希望在这个实验中清楚地听到振铃，以便比较线性相位和最小相位情况。

N = 8;     % filter order
Rp = 0.5;  % passband ripple (dB)
Rs = 60;   % stopband ripple (-dB)
Fs = 8192; % default sampling rate (Windows Matlab)
Fp = 2000; % passband end
Fc = 2200; % stopband begins [gives order 8]
[B,A] = ellip(N, Rp, Rs, Fp/(0.5*Fs)); % Matlab
fvtool(B,A);

1.2 构建两种滤波器

将上述椭圆滤波器(B,A)应用在冲击信号，得到有限冲击响应$h1$。对$h1$再应用一次椭圆滤波器(B,A)得到最小相位滤波器的响应$hmp$。将$h1$翻转后，再应用一次椭圆滤波器(B,A)得到零相位滤波器的响应$hzp$。因为$hzp$已经是对称的，不需要再翻转一次。通过这个示例，可以清晰地听到零相位滤波器的振铃效应。

Ns = 4096; % number of samples in impulse responses
imp = [1,zeros(1,Ns/2-1)]; % or 'h1=impz(B,A,Ns/2-1)'
h1 = filter(B,A,imp); % min-phase impulse response
% Minimum phase case:
hmp = filter(B,A,[h1,zeros(1,Ns/2)]); % apply twice
% Zero phase case:
h1r = fliplr(h1); % maximum-phase impulse response
hzp = filter(B,A,[h1r,zeros(1,Ns/2)]); % min*max=zp
% Let's hear them!
while(1)
  sound(hmp,Fs);
  pause(0.5);
  sound(hzp,Fs);
  pause(0.5);
end

2. 最小相位滤波器

2.1 最小相位的定义

对于一个稳定的线性时不变系统，如果它的所有零点都在单位圆内，则该系统被称为最小相位系统。

2.2 最小相位意味最快的衰减

与具有相同幅度谱的所有因果信号相比，最小相位信号最大程度地集中于时间0。由于这种特性，最小相位信号有时被称为最小延迟信号。

2.3 最大相位滤波器

对于一个稳定的线性时不变系统，如果它的所有零点都在单位圆外，则该系统被称为最大相位。最大相位的例子就是上一篇博客中提到的纯延迟和全通滤波器。

1. 纯延迟，传递函数为$H(z)=z^{-n}$，其所有极点都在$z=0$处，而零点都在$z=\infty$
2. 全通滤波器，传递函数为
   $$H(z)=\frac{c+z^{-1}}{1+c{z}^{-1}}$$ $$c=\frac{tan(\pi \ast fc/fs)-1}{tan(\pi \ast fc/fs)+1}$$其中fc为截止频率，fs为采样频率，并且$|c|<1$。因此所有的极点都在单位圆内，而所有零点都在单位圆外，是最大相位滤波器。

2.4 最小相位/全通分解

每一个单位圆上没有零点的因果稳定滤波器$H(z)$都可以被分解成一个因果稳定全通滤波器和最小相位滤波器的级联:
$$H(z)=H_{hmp}(z)H_{ap}(z)$$

举例，传递函数$H(z)$拥有一个在单位圆外的零点$a,|a|>1$，$H(z)=H_1(z)(1-a{z}^{-1})$。可以将这个传递函数按以下方式分解，其中$H_2(z)$是最小相位滤波器起 $S_1(z)$是全通滤波器
$$H(z)=H_2(z)S_1(z)$$
$$H(z)=[H_1(z)(1-a^{-1}{z}^{-1})][\frac{(1-a{z}^{-1})}{(1-a^{-1}{z}^{-1})}]$$

因此，我们将$H(z)$分解为$H_2(z)$和$S_1(z)$的乘积，其中最大相位的零点转移到单位圆内而成为最小相位(从$z=a$到$z=1/a$)。稳定的全通滤波器$S_1(z)$，该滤波器由原始最大相位的零点$a$和新极点$z=1/a$组成(它抵消了$H_2(z)$)在$z=1/a$处的零点。现在可以对$H(z)$中的每个最大相位零点重复此过程。

总之，我们可以在不改变幅频响应的情况下，将最大相位零从传递函数中提取出来，并用最小相位零替换它们。这种修改相当于将一个稳定的全通滤波器与原始滤波器串联在一起，其中全通滤波器抵消了最大相位零点并引入了最小相位零点。

3. 设计最小相位系统

最小相位滤波器设计通常需要创建最小相位所需频率响应$H(e^{j\omega })$从给定的幅度响应$|H(e^{j\omega })|$。从上述2.4可知，任何滤波器传递函数$H(z)$都可以实现最小相位。然而，对这么大的多项式进行因式分解是不切实际的。近似的“非参数”方法基于复倒谱的性质，即频谱中的每个最小相位零都会产生倒谱中的因果指数，而每个最大相位零对应于倒谱中的反因果指数。因此，通过计算倒谱并将反因果指数转换为因果指数，相应的频谱被非参数转换为最小相位形式。

matlab函数mps.m实现了这种方法。它对平滑的期望频率响应曲线效果最好，但原则上，通过简单地扩大所使用的FFT大小，可以使误差任意小。具体来说，对数量级频率响应的逆FFT不应该在时域中“缠绕”(可以忽略不计的“时间混叠”)。

在使用“相敏”滤波器设计软件(如matlab中的invfreqz)基于幅度频率响应规范设计数字滤波器时，使用mps等工具是很重要的。换句话说，当要求相敏滤波器设计软件设计因果、稳定的零相位滤波器时，通常会得到较差的结果。作为一般规则，当相位无关紧要时，要求最小相位。

一个实用的注意事项是，不稳定的递归滤波器设计通常可以通过简单地给期望的脉冲响应增加更多的延迟来稳定(在期望的相位响应上增加一个负倾斜的线性相位)。例如，Matlab中的Steiglitz-McBride算法(stmcb)是一个相位敏感的IIR滤波器设计函数，它接受期望的脉冲响应，而Matlab的invfreqz(可以选择性地迭代到Steiglitz-McBride解决方案)接受复杂的期望频率响应。

4. 文献参考：

https://www.dsprelated.com/freebooks/filters/Minimum_Phase_Filters.html