中值滤波、低通与高通滤波

1、 中值滤波

首先给出结论，中值滤波，例如说设置窗长为5个点的均值滤波，属于低通滤波。这点很容易理解，假设窗长为无限长，原始信号就变为了直流分量，频率为0。因此，均值滤波属于低通滤波，中值滤波也是一样的道理，也属于低通滤波。

2、低通滤波

我们接下来细细探究为何均值滤波属于低通滤波？
首先，例如我们得到一段随机信号，这里我们用matlab生成。

close all
clear
clc
Fs=1000;
T=1/Fs;
L=1000;
t=(0:L-1)*T;
x=0.7*sin(2*pi*20*t)+sin(2*pi*120*t);%原始信号由20Hz和120Hz的2段信号组成
% x=x+2*randn(size(t));%可添加随机白噪声，也可不添加
figure(1)
plot(t,x)
hold on

原始信号图片如下：

低通滤波

% 4阶低通滤波
[b,a]=butter(4,50/(Fs/2),'low');
y1=filter(b,a,x);
plot(t,y1,'LineWidth',2)

可以看到，120Hz的信号成分基本完全被滤除了。

3、均值滤波

%均值滤波
y2=smooth(x,5);
figure(2)
plot(t,x)
hold on
plot(t,y2,'LineWidth',2)

可以看到，均值滤波并没有将120Hz的成分所滤除。

3.1 低通滤波数学表达

这里我们看到了，低通滤波用的2行代码是这样的

% 4阶低通滤波
[b,a]=butter(4,50/(Fs/2),'low');
y1=filter(b,a,x);

这两行代码对原数组进行了何种操作呢？
首先，butter是用来生成零极点的滤波器系数。这里我们说明一下，滤波器通常分为FIR有限冲击响应滤波器和IIR无限冲击响应滤波器。
IIR滤波器的转移函数为
$$H(z)=\frac{{\sum }_{r=0}^M{b}_r{z}^{-r}}{1+{\sum }_{k=1}^N{a}_k{z}^{-k}}$$
FIR滤波器的转移函数为
$$H(z)=\sum _{n=0}^{N-1}h(n)z^{-n}$$
这里我们能看出来，这里使用的带通滤波器属于IIR滤波器，中值滤波属于FIR滤波。
在matlab中可以观察到

a=[1,-3.1806,3.8612,-2.1122,0.4383];
b=[-0.0004,0.0017,0.0025,0.0017,0.0004];

在filter函数中做了何种处理呢？我们来看下官方matlab对filter.m的解释

%FILTER One-dimensional digital filter.
%   Y = FILTER(B,A,X) filters the data in vector X with the
%   filter described by vectors A and B to create the filtered
%   data Y.  The filter is a "Direct Form II Transposed"
%   implementation of the standard difference equation:
%
%   a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
%                         - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
%
%   If a(1) is not equal to 1, FILTER normalizes the filter
%   coefficients by a(1). 

这里我们需要理解注视中间的公式是如何得出来的。
我们知道，在Z变换中
$$H(z)=\frac{Y(Z)}{X(Z)}$$
因此，由IIR的转移函数，我们可知
$$Y(z)[1+\sum _{k=1}^{N}a(k)z^{-k}]=X(z)[b(0)+\sum _{r=1}^{M}b(r)z^{-r}]$$
进一步可知
$$Y(z)=-Y(z)\sum _{k=1}^{N}a(k)z^{-k}+X(z)\sum _{r=0}^{M}b(r)z^{-r}$$
由Z的逆变换，最后可推出
$$y(n)=-\sum _{k=1}^{N}a(k)y(n-k)+\sum _{r=0}^{M}b(r)x(n-r)$$
因此，在filter中，进行的操作如下

%这里需要注意在matlab中数组下标由1开始，而前面公式推导中默认下标从0开始
y(1)=x(1);
y(2)=-(a(2)*y(1))+b(1)*x(2)+b(2)*x(1)=0.000321;
y(3)=-(a(2)*y(2)+a(3)*y(1))+b(1)*x(3)+b(2)*x(2)+b(3)*x(1)=0.0028;
y(4)=-(a(2)*y(3)+a(3)*y(2)+a(4)*y(1))+b(1)*x(4)+b(2)*x(3)+b(3)*x(2)+b(4)*x(1)=0.0120
......

我们将由b,a组成的滤波器系统的频率响应图画出

freqz(b,a,[],Fs)

可看出，50Hz以上的频率成分被衰减了。
这里，我们再看下中值滤波的形式。
我们讨论两种形式，一种只利用前面时刻的数据进行中值滤波，一种还利用未来时刻的数据进行滤波。
$$y(n)=\frac{x(n)+x(n-1)+x(n-2)+x(n-3)+x(n-4)}{5};$$
由Z变换得到
$$Y(z)=\frac{X(z)+X(z)z^{-1}+X(z)z^{-2}+X(z)z^{-3}+X(z)z^{-4}}{5}$$
进一步得到
$$H(z)=0.2+0.2\ast z^{-1}+0.2\ast z^{-2}+0.2\ast z^{-3}+0.2\ast z^{-4}$$
我们可以得出均值滤波的系数为

b=[0.2,0.2,0.2,0.2,0.2];
a=[1,0,0,0,0];

我们同样绘制出频率响应曲线图，我们可以另开一个小标题专门讲下。

如何根据传递函数绘制幅频响应曲线图

我们现在再次回到传递函数，
$$H(z)=\frac{{\sum }_{r=0}^M{b}_r{z}^{-r}}{1+{\sum }_{k=1}^N{a}_k{z}^{-k}}$$
对分子和分母分别做因式分解，得到
$$H(z)=g{z}^{N-M}\frac{{\prod }_{r=1}^M(z-z_r)}{{\prod }_{k=1}^N(z-p_k)}$$
g为系统的增益因子，g=b(0)。pk为极点，zr为系统的零点。求出零极点后可根据零极点图做出系统的幅频响应曲线。
令$z=e^{jw}$，即令z在单位圆上取值，可进行变换
$$H(e^{jw})=g{e}^{j(N-M)w}\frac{{\prod }_{r=1}^M(e^{jw}-z_r)}{{\prod }_{k=1}^N(e^{jw}-p_k)}$$
式中，$e^{jw}$是单位圆上的某一个点，$e^{jw}$可以看作是原点到$e^{jw}$的向量，$z_r$是平面上的零点，$e^{jw}-z_r$可看作是由零点$z_r$到$e^{jw}$的向量。分母中的极点情况类似，因此，我们可以得到系统幅频响应的几何解释
$$|H(e^{jw})|=g\frac{{\prod }_{r=1}^M|e^{jw}-z_r|}{{\prod }_{k=1}^N|e^{jw}-p_k|}$$
省略的中间的某一项模值为0.
这里懒得打公式了，直接在ipad上作图说明一下。
首先为了简化，我们讲5点均值滤波改为2点均值滤波。

纠正下上面图片的错误，我们真实取的周期根据对称性应该是$[-\pi ,\pi ]$，因此均值滤波器还是低通滤波器。
===分割线
这里在提供一种思路
例如我们考虑一个简单的例子，在输入值上取2点的平均：
$$y[n]=\frac{1}{2}(x[n]+x[n-1])$$
在这种情况下，
$$h[n]=\frac{1}{2}(\delta [n]+\delta [n-1])$$
根据我们已知的z变换的定义，$H(z)={\sum }_{n=-\infty }^{\infty }h[n]z^{-n}$,将$z=e^{jw}$代入推导出
$$H(e^{jw})=\sum _{n=-\infty }^{\infty }h[n]e^{-jwn}$$
根据上述公式，代入可推导出2点均值滤波的系统频率响应为
$$\begin{gathered} H(e^{jw})=h[0]+h[1]e^{-jw} \\ =\frac{1}{2}(1+e^{-jw})=\frac{1}{2}{e}^{-jw/2}(e^{jw/2+e^{-jw/2}}) \\ =e^{-jw/2}cos(w/2) \end{gathered}$$
因此，可以知道幅频响应曲线如下图所示。

从上图中我们可以看出来，均值滤波严格来说算是属于低通滤波器的一种，但是对于高频成分也没有完全滤除。实际上，均值滤波属于频率成形滤波器，而我们常说的高低通滤波则属于频率选择滤波器。
定义：
频率成形滤波器：用于改变频谱形状的线性时不变系统往往成为频率成形滤波器。
频率选择滤波器：专门设计成基本上无失真低通过某些频率，而显著衰减或消除另一些频率的系统成为频率选择性滤波器。