实时零相位滤波的神话（1）

做控制的人大概都梦想做到输出和输入信号保持完全同步，相移为0。如果能做到，那该多酷。

MATLAB有个神奇的函数filtfilt，可以对数据做离线的滤波，实现零相移。原理就是先做一个方向的滤波，比如先forward 滤波，然后把滤波后的序列逆序，再用同一个滤波器做backward滤波，滤波得到的序列最后再逆序，得到最终结果。尝试了一下，滤波效果那真是赏心悦目。可惜，不能实现在线信号的实时零相移滤波。怎么做呢？焦虑之际，在百度上搜到这个文章：http://blog.csdn.net/shenziheng1/article/details/53415642，作者的手段让我万分佩服：既然相移导致的group delay不为0，那咱就设计个group delay与其完全相反的全通滤波器，二者级联，groupdelay两两相消，不就搞定了吗？多么完美！！

感谢这位博主，在相关文章里把代码给的相当完整，某些略去的地方其实比较容易补全。按照这个思路做下去，确实得到了相移接近0 的结果，图片如下：

但是当我兴冲冲地把low pass 和all pass级联进行滤波，得到了让人郁闷的玩意：滤波结果是震荡的，数据越来越大。。。。痛苦地找了好久，某篇英文文章说，级联allpass滤波器，只能增大相移，尤其是在将iir整成线性相位滤波器时的典型应用时，group delay在通带内虽然可以保持恒定，但只能让group delay变大，即输出结果比输入delay更多个samples。这个震荡，应该是all pass的零极点与单位圆的关系导致系统不稳定造成的。

iir+all pass常用来进行相位均衡，这个均衡并不是要做到零相位，而是使得通带内的信号保持固定的group delay。下面是从某个大学课件上搞来的iir+all pass的代码，讲解了相位均衡的常规应用：

t = 0:0.001:1;


w = pi*t; % normalized frequencies
alpha = 0.5; % LPF parameter
blp = (1-alpha)/2*[1 1]; % numerator coefficients
alp = [1 -alpha]; % denominator coefficients
hlp = freqz(blp,alp,w); % compute DTFT
glp = grpdelay(blp,alp,w); % compute group delay


N = 4; % all-pass filter order
F = w(1:501)/pi; % normalized frequencies
edges = [0 1/2]; % band-edge frequencies
Gd = max(glp)-glp(1:501); % desired group-delays of APF (>0)
[bap,aap] = iirgrpdelay(N,F,edges,Gd); % make all-pass filter
hap = freqz(bap,aap,w); % compute DTFT
gap = grpdelay(bap,aap,w); % compute group-delay

b = conv(blp,bap); % product of numerators
a = conv(alp,aap); % product of denominators
h = freqz(b,a,w); % compute DTFT
g = grpdelay(b,a,w); % compute group delay
subplot(3,1,1)
plot(w,abs(h));
ylabel('magnitude response');
subplot(3,1,2);
plot(w,unwrap(angle(h)));
ylabel('phase response (rad)');
subplot(3,1,3);
plot(w,g);
xlabel('normalized freq (rad/sample)');
ylabel('group delay (samples)');

已经确定了，all pass 级联不能让group delay 变成0，反而是更大。于是，我的第一次充满勇气的寻找零相位实时滤波器行动宣告破产。