维纳（Wiener）滤波及Matlab代码

文章目录

维纳（Wiener）滤波
模型结构
使用条件
原理公式推导
仿真分析——Matlab代码
- 一、参考信号d(n)d(n)d(n)为原始信号s(n)s(n)s(n)
- 二、参考信号d(n)d(n)d(n)为加性高斯白噪声v(n)v(n)v(n)
- 三、参考信号d(n)d(n)d(n)为输入信号自身x(n)x(n)x(n)
总结

维纳（Wiener）滤波

Wiener滤波的核心目标就是使均方误差最小，从而可以推导出维纳-霍夫方程。

在信号处理中，滤波器可以分为FIR和IIR两种，在这里主要介绍维纳FIR滤波器，如下图1所示。换句话说，就是要想方设法求出最优滤波器的系数。

图1 Wiener滤波的横向滤波器

模型结构

首先，先来看一下wiener滤波器的一般结构，如下图2所示。

图2 wiener滤波的原理框图

其中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[e2(n)]J(h)=E[{e^2}(n)]J(h)=E[e2(n)]

为了使均方误差最小，需要对JJJ进行求导，并让导数为0，即可得到最佳滤波器的系数了。

∇hJ=−2E[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))]=0E\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)rx(i−k)=rxd(−k)\sum_{i=-\infty}^{\infty} h(i)r_{x}(i-k)=r_{x d}(-k) i=−∞∑∞h(i)rx(i−k)=rxd(−k)

针对MMM阶FIR滤波器，维纳-霍夫方程（Wiener-Hopf）为∑i=0M−1h(i)rx(i−k)=rxd(−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)，那么，写成矩阵的形式就是
Rh=rxdh=R−1rxd\begin{array}{l} Rh={r}_{x d} \\ h=R^{-1}{r}_{x d} \end{array} Rh=rxdh=R−1rxd
其中RRR是信号x(n)x(n)x(n)的自相关矩阵，rxdr_{x d}rxd是信号x(n)和d(n)x(n)和d(n)x(n)和d(n)的互相关矩阵。

仿真分析——Matlab代码

由于Wiener滤波器的参数求解过程中必须要用到参考信号，所以这里仿真采用的信号Signal为
s=A∗sin(2πf1t)+B∗sin(2πf2t)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,f1=1000,f2=2000A = 10,B = 15,{f_1} = 1000,{f_2} = 2000A=10,B=15,f1=1000,f2=2000 。
根据采样定理，这里假定采样频率fs=100000{f_s} = 100000fs=100000，采样间隔T=1/fsT = 1/{f_s}T=1/fs，则sa(t)∣t=nT=sa(nT)(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)。
对于不同的nnn值，s(n)s(n)s(n)是一个有序的数字序列:s(n)={sa(0),sa(T),sa(2T),⋯}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)也可以称作参考信号。有以下几种不同的情况：

参考信号d(n)d(n)d(n)为原始信号s(n)s(n)s(n)
参考信号d(n)d(n)d(n)为加性高斯白噪声v(n)v(n)v(n)
参考信号d(n)d(n)d(n)为输入信号自身x(n)x(n)x(n)

下面对三种不同的情形进行仿真与对比分析, 其中滤波器系数长度N均为100。
（全部三种完整的MATLAB代码整合版在我的资源“维纳（Wiener）滤波及Matlab代码”中）

一、参考信号d(n)d(n)d(n)为原始信号s(n)s(n)s(n)

%% 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:Nfor j=1:NmRxx(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:NmRxd(i)=Rxd(N-1+i); % 1*N维度
endh = 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(n)d(n)d(n)为加性高斯白噪声v(n)v(n)v(n)

只要将代码的开头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(n)d(n)d(n)为输入信号自身x(n)x(n)x(n)

只要将代码的开头d换成x即可。

滤波后信噪比SNR_after = -0.812dB；滤波后较滤波前信噪比提高了-3.812dB。

总结

参考信号选为原始信号时的滤波效果最好。虽然在第三种方案中，参考信号选为自身信号时的滤波效果不好，第三种Wiener滤波器模型对于许多应用仍然具有很大的实用价值, 因为在许多实际应用中, 既无法提前获知系统的期望响应, 也不具备获得与输入信号高相关噪声的条件。