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

中值滤波

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

低通滤波

我们接下来细细探究为何均值滤波属于低通滤波?
首先,例如我们得到一段随机信号,这里我们用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

原始信号图片如下:
中值滤波、低通与高通滤波_第1张图片
低通滤波

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

中值滤波、低通与高通滤波_第2张图片
可以看到,120Hz的信号成分基本完全被滤除了。

均值滤波

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

中值滤波、低通与高通滤波_第3张图片
可以看到,均值滤波并没有将120Hz的成分所滤除。

低通滤波数学表达

这里我们看到了,低通滤波用的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)

中值滤波、低通与高通滤波_第4张图片
可看出,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)=5(x(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)=5(X(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];

我们同样绘制出频率响应曲线图

你可能感兴趣的