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 ) = ∑ r = 0 M b r z − r 1 + ∑ k = 1 N a k z − k H(z)=\frac{\sum_{r=0}^{M}b_{r}z^{-r}}{1+\sum_{k=1}^{N}a_{k}z^{-k}} H(z)=1+k=1Nakzkr=0Mbrzr
FIR滤波器的转移函数为
H ( z ) = ∑ n = 0 N − 1 h ( n ) z − n H(z)=\sum_{n=0}^{N-1}h(n)z^{-n} H(z)=n=0N1h(n)zn
这里我们能看出来,这里使用的带通滤波器属于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 ) = Y ( Z ) X ( Z ) H(z)=\frac{Y(Z)}{X(Z)} H(z)=X(Z)Y(Z)
因此,由IIR的转移函数,我们可知
Y ( z ) [ 1 + ∑ k = 1 N a ( k ) z − k ] = X ( z ) [ b ( 0 ) + ∑ r = 1 M b ( r ) z − r ] 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)[1+k=1Na(k)zk]=X(z)[b(0)+r=1Mb(r)zr]
进一步可知
Y ( z ) = − Y ( z ) ∑ k = 1 N a ( k ) z − k + X ( z ) ∑ r = 0 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} Y(z)=Y(z)k=1Na(k)zk+X(z)r=0Mb(r)zr
由Z的逆变换,最后可推出
y ( n ) = − ∑ k = 1 N a ( k ) y ( n − k ) + ∑ r = 0 M b ( r ) x ( n − r ) y(n)=-\sum_{k=1}^{N}a(k)y(n-k)+\sum_{r=0}^{M}b(r)x(n-r) y(n)=k=1Na(k)y(nk)+r=0Mb(r)x(nr)
因此,在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 ) = x ( n ) + x ( n − 1 ) + x ( n − 2 ) + x ( n − 3 ) + x ( n − 4 ) 5 ; y(n)=\frac{x(n)+x(n-1)+x(n-2)+x(n-3)+x(n-4)}{5}; y(n)=5x(n)+x(n1)+x(n2)+x(n3)+x(n4);
由Z变换得到
Y ( z ) = X ( z ) + X ( z ) z − 1 + X ( z ) z − 2 + X ( z ) z − 3 + X ( z ) z − 4 5 Y(z)=\frac{X(z)+X(z)z^{-1}+X(z)z^{-2}+X(z)z^{-3}+X(z)z^{-4}}{5} Y(z)=5X(z)+X(z)z1+X(z)z2+X(z)z3+X(z)z4
进一步得到
H ( z ) = 0.2 + 0.2 ∗ z − 1 + 0.2 ∗ z − 2 + 0.2 ∗ z − 3 + 0.2 ∗ z − 4 H(z)=0.2+0.2*z^{-1}+0.2*z^{-2}+0.2*z^{-3}+0.2*z^{-4} H(z)=0.2+0.2z1+0.2z2+0.2z3+0.2z4
我们可以得出均值滤波的系数为

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

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

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

我们现在再次回到传递函数,
H ( z ) = ∑ r = 0 M b r z − r 1 + ∑ k = 1 N a k z − k H(z)=\frac{\sum_{r=0}^{M}b_{r}z^{-r}}{1+\sum_{k=1}^{N}a_{k}z^{-k}} H(z)=1+k=1Nakzkr=0Mbrzr
对分子和分母分别做因式分解,得到
H ( z ) = g z N − M ∏ r = 1 M ( z − z r ) ∏ k = 1 N ( z − p k ) H(z)=gz^{N-M}\frac{\prod_{r=1}^{M}(z-z_{r})}{\prod_{k=1}^{N}(z-p_{k})} H(z)=gzNMk=1N(zpk)r=1M(zzr)
g为系统的增益因子,g=b(0)pk为极点,zr为系统的零点。求出零极点后可根据零极点图做出系统的幅频响应曲线。
z = e j w z=e^{jw} z=ejw,即令z在单位圆上取值,可进行变换
H ( e j w ) = g e j ( N − M ) w ∏ r = 1 M ( e j w − z r ) ∏ k = 1 N ( e j w − p k ) H(e^{jw})=ge^{j(N-M)w}\frac{\prod_{r=1}^{M}(e^{jw}-z_{r})}{\prod_{k=1}^{N}(e^{jw}-p_{k})} H(ejw)=gej(NM)wk=1N(ejwpk)r=1M(ejwzr)
式中, e j w e^{jw} ejw是单位圆上的某一个点, e j w e^{jw} ejw可以看作是原点到 e j w e^{jw} ejw的向量, z r z_{r} zr是平面上的零点, e j w − z r e^{jw}-z_{r} ejwzr可看作是由零点 z r z_{r} zr e j w e^{jw} ejw的向量。分母中的极点情况类似,因此,我们可以得到系统幅频响应的几何解释
∣ H ( e j w ) ∣ = g ∏ r = 1 M ∣ e j w − z r ∣ ∏ k = 1 N ∣ e j w − p k ∣ |H(e^{jw})|=g\frac{\prod_{r=1}^{M}|e^{jw}-z_{r}|}{\prod_{k=1}^{N}|e^{jw}-p_{k}|} H(ejw)=gk=1Nejwpkr=1Mejwzr
省略的中间的某一项模值为0.
这里懒得打公式了,直接在ipad上作图说明一下。
首先为了简化,我们讲5点均值滤波改为2点均值滤波。
在这里插入图片描述
纠正下上面图片的错误,我们真实取的周期根据对称性应该是 [ − π , π ] [-\pi,\pi] [π,π],因此均值滤波器还是低通滤波器。
===分割线
这里在提供一种思路
例如我们考虑一个简单的例子,在输入值上取2点的平均:
y [ n ] = 1 2 ( x [ n ] + x [ n − 1 ] ) y[n]=\frac{1}{2}(x[n]+x[n-1]) y[n]=21(x[n]+x[n1])
在这种情况下,
h [ n ] = 1 2 ( δ [ n ] + δ [ n − 1 ] ) h[n]=\frac{1}{2}(\delta[n]+\delta[n-1]) h[n]=21(δ[n]+δ[n1])
根据我们已知的z变换的定义, H ( z ) = ∑ n = − ∞ ∞ h [ n ] z − n H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n} H(z)=n=h[n]zn,将 z = e j w z=e^{jw} z=ejw代入推导出
H ( e j w ) = ∑ n = − ∞ ∞ h [ n ] e − j w n H(e^{jw})=\sum_{n=-\infty}^{\infty}h[n]e^{-jwn} H(ejw)=n=h[n]ejwn
根据上述公式,代入可推导出2点均值滤波的系统频率响应为
H ( e j w ) = h [ 0 ] + h [ 1 ] e − j w = 1 2 ( 1 + e − j w ) = 1 2 e − j w / 2 ( e j w / 2 + e − j w / 2 ) = e − j w / 2 c o s ( w / 2 ) 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) H(ejw)=h[0]+h[1]ejw=21(1+ejw)=21ejw/2(ejw/2+ejw/2)=ejw/2cos(w/2)
因此,可以知道幅频响应曲线如下图所示。

在这里插入图片描述

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

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐