中值滤波、低通与高通滤波
首先给出结论,中值滤波,例如说设置窗长为5个点的均值滤波,属于低通滤波。这点很容易理解,假设窗长为无限长,原始信号就变为了直流分量,频率为0。因此,均值滤波属于低通滤波,中值滤波也是一样的道理。
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=1Nakz−k∑r=0Mbrz−r
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=0∑N−1h(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
)
=
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=1∑Na(k)z−k]=X(z)[b(0)+r=1∑Mb(r)z−r]
进一步可知
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=1∑Na(k)z−k+X(z)r=0∑Mb(r)z−r
由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=1∑Na(k)y(n−k)+r=0∑Mb(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
)
=
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(n−1)+x(n−2)+x(n−3)+x(n−4);
由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)z−1+X(z)z−2+X(z)z−3+X(z)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.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
我们可以得出均值滤波的系数为
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=1Nakz−k∑r=0Mbrz−r
对分子和分母分别做因式分解,得到
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)=gzN−M∏k=1N(z−pk)∏r=1M(z−zr)
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(N−M)w∏k=1N(ejw−pk)∏r=1M(ejw−zr)
式中,
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}
ejw−zr可看作是由零点
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)∣=g∏k=1N∣ejw−pk∣∏r=1M∣ejw−zr∣
省略的中间的某一项模值为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[n−1])
在这种情况下,
h
[
n
]
=
1
2
(
δ
[
n
]
+
δ
[
n
−
1
]
)
h[n]=\frac{1}{2}(\delta[n]+\delta[n-1])
h[n]=21(δ[n]+δ[n−1])
根据我们已知的z变换的定义,
H
(
z
)
=
∑
n
=
−
∞
∞
h
[
n
]
z
−
n
H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}
H(z)=∑n=−∞∞h[n]z−n,将
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]e−jwn
根据上述公式,代入可推导出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]e−jw=21(1+e−jw)=21e−jw/2(ejw/2+e−jw/2)=e−jw/2cos(w/2)
因此,可以知道幅频响应曲线如下图所示。
从上图中我们可以看出来,均值滤波严格来说算是属于低通滤波器的一种,但是对于高频成分也没有完全滤除。实际上,均值滤波属于频率成形滤波器,而我们常说的高低通滤波则属于频率选择滤波器。
定义:
频率成形滤波器:用于改变频谱形状的线性时不变系统往往成为频率成形滤波器。
频率选择滤波器:专门设计成基本上无失真低通过某些频率,而显著衰减或消除另一些频率的系统成为频率选择性滤波器。
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)