赞
踩
Wiener滤波的核心目标就是使均方误差最小,从而可以推导出维纳-霍夫方程。
在信号处理中,滤波器可以分为FIR和IIR两种,在这里主要介绍维纳FIR滤波器,如下图1所示。换句话说,就是要想方设法求出最优滤波器的系数。
首先,先来看一下wiener滤波器的一般结构,如下图2所示。
其中 s ( n ) s(n) s(n)是输入的原始信号,经过噪声信道 v ( n ) v(n) v(n)后,变成了 x ( n ) x(n) x(n),滤波器后输出得到 s ′ ( n ) s'(n) s′(n),期望响应是 d ( n ) d(n) d(n) ,那误差 e ( n ) = d ( n ) − s ′ ( n ) e(n)=d(n)-s'(n) e(n)=d(n)−s′(n) ,滤波器的系数为 h ( n ) h(n) h(n)。
对于信号 s ( n ) s(n) s(n)和噪声 v ( n ) v(n) v(n)的混合体 x ( n ) = s ( n ) + v ( n ) x(n)=s(n)+v(n) x(n)=s(n)+v(n),按照均方误差最小的准则,从 x ( t ) x(t) x(t)中分离出信号 s ( t ) s(t) s(t)的理论,称为维纳滤波理论。
这里要着重说明的一点是,期望响应 d ( n ) d(n) d(n)一般上来说是对原始信号的 s ( n ) s(n) s(n)估计,后面会讲到别的情况。如果你想要对一个未知的信号进行滤波,用wiener滤波的方法是不行的。因为,在设计滤波器系数的时候,必须用到期望信号 d ( n ) d(n) d(n)。所以,必须要知道 d ( n ) d(n) d(n),或者 d ( n ) d(n) d(n)的一些其他的特征。
1、输入信号是宽平稳信号。那么宽平稳是什么呢?通俗的来说就是与时间无关的信号,或者与时间的起点无关,只与时间间隔有关。
2、必须知道期望信号,或者它的一些信号特征才行(具体看下面的公式推导部分)。
为了简便运算,假设所使用的信号是实信号。
输出信号: s ′ ( n ) = h ( n ) ∗ x ( n ) = ∑ k = − ∞ + ∞ h ( n ) x ( n − k ) s'(n)=h(n)*x(n)=\sum_{k=-\infty}^{+\infty} h(n) x(n-k) s′(n)=h(n)∗x(n)=∑k=−∞+∞h(n)x(n−k)
误差: e ( n ) = d ( n ) − s ′ ( n ) = = d ( n ) − ∑ k = − ∞ + ∞ h ( n ) x ( n − k ) e(n)=d(n)-s'(n)==d(n)-\sum_{k=-\infty}^{+\infty} h(n)x(n-k) e(n)=d(n)−s′(n)==d(n)−∑k=−∞+∞h(n)x(n−k)
均方误差: J ( h ) = E [ e 2 ( n ) ] J(h)=E[{e^2}(n)] J(h)=E[e2(n)]
为了使均方误差最小,需要对 J J J进行求导,并让导数为0,即可得到最佳滤波器的系数了。
∇ h J = − 2 E [ x ( n − k ) e ( n ) ] {\nabla _h}J = - 2E[ x(n - k)e(n)] ∇hJ=−2E[x(n−k)e(n)]
所以,
E
[
x
(
n
−
k
)
e
(
n
)
]
=
0
(
1
)
E[x(n - k)e(n)]=0\quad\quad\quad\quad\quad(1)
E[x(n−k)e(n)]=0(1)
E [ x ( n − k ) ( d ( n ) − ∑ i = − ∞ ∞ h ( i ) x ( n − i ) ) ] = 0 E\left[x(n-k)\left(d(n)-\sum_{i=-\infty}^{\infty} h(i)x(n-i)\right)\right]=0 E[x(n−k)(d(n)−i=−∞∑∞h(i)x(n−i))]=0
∑ i = − ∞ ∞ h ( i ) E [ x ( n − k ) x ( n − i ) ] = E [ x ( n − k ) d ( n ) ] \sum_{i=-\infty}^{\infty} h(i) E\left[x(n-k) x(n-i)\right]=E\left[x(n-k)d(n)\right] i=−∞∑∞h(i)E[x(n−k)x(n−i)]=E[x(n−k)d(n)]
∑ i = − ∞ ∞ h ( i ) r x ( i − k ) = r x d ( − k ) \sum_{i=-\infty}^{\infty} h(i)r_{x}(i-k)=r_{x d}(-k) i=−∞∑∞h(i)rx(i−k)=rxd(−k)
针对
M
M
M阶FIR滤波器,维纳-霍夫方程(Wiener-Hopf)为
∑
i
=
0
M
−
1
h
(
i
)
r
x
(
i
−
k
)
=
r
x
d
(
−
k
)
\sum_{i=0}^{M-1} h(i)r_{x}(i-k)=r_{x d}(-k)
∑i=0M−1h(i)rx(i−k)=rxd(−k),那么,写成矩阵的形式就是
R
h
=
r
x
d
h
=
R
−
1
r
x
d
其中
R
R
R是信号
x
(
n
)
x(n)
x(n)的自相关矩阵,
r
x
d
r_{x d}
rxd是信号
x
(
n
)
和
d
(
n
)
x(n)和d(n)
x(n)和d(n)的互相关矩阵。
由于Wiener滤波器的参数求解过程中必须要用到参考信号,所以这里仿真采用的信号Signal为
s
=
A
∗
s
i
n
(
2
π
f
1
t
)
+
B
∗
s
i
n
(
2
π
f
2
t
)
s = A*sin\left( {2\pi {f_1}t} \right){\rm{ }} + {\rm{ }}B*sin\left( {2\pi {f_2}t} \right)
s=A∗sin(2πf1t)+B∗sin(2πf2t) ,
即为两个叠加的正弦信号。其中
A
=
10
,
B
=
15
,
f
1
=
1000
,
f
2
=
2000
A = 10,B = 15,{f_1} = 1000,{f_2} = 2000
A=10,B=15,f1=1000,f2=2000 。
根据采样定理,这里假定采样频率
f
s
=
100000
{f_s} = 100000
fs=100000,采样间隔
T
=
1
/
f
s
T = 1/{f_s}
T=1/fs,则
s
a
(
t
)
∣
t
=
n
T
=
s
a
(
n
T
)
(
0
≤
n
≤
999
)
{s_a}(t)\left| {_{t = nT}} \right. = {s_a}(nT){\rm{ }}(0 \le n \le 999)
sa(t)∣t=nT=sa(nT)(0≤n≤999)。
对于不同的
n
n
n值,
s
(
n
)
s(n)
s(n)是一个有序的数字序列:
s
(
n
)
=
{
s
a
(
0
)
,
s
a
(
T
)
,
s
a
(
2
T
)
,
⋯
}
s(n) = \left\{ {{s_a}(0),{s_a}(T),{s_a}(2T), \cdots } \right\}
s(n)={sa(0),sa(T),sa(2T),⋯}。信号传输过程中,信道中的噪声为加性高斯白噪声。原始信噪比定义为SNR=3。
上面提到期望信号 d ( n ) d(n) d(n)是必不可少的,所以,对期望信号的设定也会决定结果的好坏。所以, d ( n ) d(n) d(n)也可以称作参考信号。有以下几种不同的情况:
下面对三种不同的情形进行仿真与对比分析, 其中滤波器系数长度N均为100。
(全部三种完整的MATLAB代码整合版在我的资源“维纳(Wiener)滤波及Matlab代码”中)
%% wiener filter for different d(n) %% StarryHuang 2021.1.9 close all; clear; clc; %% 信号产生,对原始信号进行采样 A=10; % 信号的幅值 B=15; % 信号的幅值 f1=1000; % 信号1的频率 f2=2000; % 信号2的频率 fs=10^5; % 采样频率 t=0:999; % 采样点t = [0,999],长度1000 M = length(t); % 信号长度 s=A*sin(2*pi*f1*t/fs) + B*sin(2*pi*f2*t/fs); % 采样正弦波信号为Signal SNR = 3; % 初始信噪比 x=awgn(s,SNR,'measured'); %观测信号 x=s+v.,给正弦波信号加入信噪比为-3dB的高斯白噪声 v=x - s; % 高斯白噪声,误差信号,每次运行都不一样 %% 第一种情况——期望信号d(n)为感兴趣的原信号Signal d = s; %% 第二种情况——期望信号d(n)为Noise v % d = v; %% 第三种情况—— 期望信号d(n)为x(n-1) % d = [x(1),x(1:end-1),]; % d(n)=x(n-1) %% 维纳滤波 N = floor(length(x)*0.1); % 滤波器的阶数,向下取整 % N=100; % 维纳滤波器阶数 Rxx=xcorr(x,N-1,'biased'); % 自相关函数1*(2N-1)维度,返回一个延迟范围在[-N,N]的互相关函数序列,对称的 % 变成矩阵 N*N维度 for i=1:N for j=1:N mRxx(i,j)=Rxx(N-i+j); % N*N维度 end end %产生维纳滤波中x 方向上观测信号与期望信号d的互相关矩阵 Rxd=xcorr(x,d,N-1,'biased'); % 互相关函数1*(2N-1)维度 % 变成矩阵1*N维度 for i=1:N mRxd(i)=Rxd(N-1+i); % 1*N维度 end h = inv(mRxx)*mRxd'; % 由wiener-Hopf方程得到滤波器最优解, h是N*1维度 %% 检验wiener滤波效果 y = conv(x,h); % 滤波后的输出,长度为M+N-1,要截取前M个。 y = y(1:M); % yy = filter(h,1,x); % 用卷积或者直接用filter都可以 Py=sum(power(y,2))/M; %滤波后信号y的功率 e = d - y; % 输出减去期望等于滤波误差 J = sum(power(e,2))/M % 滤波后噪声功率 SNR_after = 10*log10((Py-J)/J) % 滤波后信噪比 db单位 imv = 10*log10((Py-J)/J/power(10,SNR/10)) % 滤波后较滤波前信噪比提高了imv dB。 %% 画图 % 原始信号x,噪声v,观测波形d figure(1), subplot(311) plot(t,s) % 输入函数 title(' Signal原信号') subplot(312) plot(t,v) % 输入函数 title(' Noise噪声') subplot(313) plot(t,x) % 输入函数 title(' X(n)观测波形') %% d = s % 期望和滤波后的信号对比 figure(2), subplot(211) plot(t, d, 'r:', t, y, 'b-','LineWidth',1); legend('期望信号','滤波后结果'); title('期望信号与滤波结果对比'); xlabel('观测点数');ylabel('信号幅度'); axis([0 1000 -50 50]) subplot(212), plot(t, e); title('输出误差'); xlabel('观测点数');ylabel('误差幅度'); axis([0 1000 -50 50]) % 滤波前后对比 figure(3), subplot(211) plot(t, x); title('维纳滤波前'); xlabel('观测点数');ylabel('信号幅度'); axis([0 1000 -50 50]) subplot(212), plot(t, y); title('维纳滤波后'); xlabel('观测点数');ylabel('误差幅度'); axis([0 1000 -50 50])
滤波后信噪比SNR_after = 13.187dB;滤波后较滤波前信噪比提高了10.4dB。
只要将代码的开头d换成v,画图函数修改为%% d=v部分即可。
% %% d = v
% figure(2)
% plot(t, d, 'r:', t, y, 'b-','LineWidth',1);
% legend('期望信号','滤波后结果'); title('期望信号与滤波结果对比');
% xlabel('观测点数');ylabel('信号幅度');
% axis([0 1000 -50 50])
%
% figure(3)
% plot(t, s, 'r:', t, x - y, 'b-','LineWidth',1);
% legend('Signal原始信号','噪声抵消后结果'); title('Signal原始信号与噪声抵消后结果对比');
% xlabel('观测点数');ylabel('信号幅度');
% axis([0 1000 -50 50])
滤波后信噪比SNR_after = 10.023dB;滤波后较滤波前信噪比提高了7dB。
只要将代码的开头d换成x即可。
滤波后信噪比SNR_after = -0.812dB;滤波后较滤波前信噪比提高了-3.812dB。
参考信号选为原始信号时的滤波效果最好。虽然在第三种方案中,参考信号选为自身信号时的滤波效果不好,第三种Wiener滤波器模型对于许多应用仍然具有很大的实用价值, 因为在许多实际应用中, 既无法提前获知系统的期望响应, 也不具备获得与输入信号高相关噪声的条件。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。