维纳（Wiener）滤波及Matlab代码

维纳（Wiener）滤波

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

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

图1 Wiener滤波的横向滤波器

模型结构

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

图2 wiener滤波的原理框图

其中$s(n)$是输入的原始信号，经过噪声信道$v(n)$后，变成了$x(n)$，滤波器后输出得到$s^{\prime }(n)$，期望响应是$d(n)$ ，那误差$e(n)=d(n)-s^{\prime }(n)$ ，滤波器的系数为$h(n)$。

对于信号$s(n)$和噪声$v(n)$的混合体$x(n)=s(n)+v(n)$，按照均方误差最小的准则，从$x(t)$中分离出信号$s(t)$的理论，称为维纳滤波理论。

这里要着重说明的一点是，期望响应$d(n)$一般上来说是对原始信号的$s(n)$估计，后面会讲到别的情况。如果你想要对一个未知的信号进行滤波，用wiener滤波的方法是不行的。因为，在设计滤波器系数的时候，必须用到期望信号$d(n)$。所以，必须要知道$d(n)$，或者$d(n)$的一些其他的特征。

使用条件

1、输入信号是宽平稳信号。那么宽平稳是什么呢？通俗的来说就是与时间无关的信号，或者与时间的起点无关，只与时间间隔有关。

2、必须知道期望信号，或者它的一些信号特征才行（具体看下面的公式推导部分）。

原理公式推导

为了简便运算，假设所使用的信号是实信号。

输出信号：$s^{\prime }(n)=h(n)\ast x(n)={\sum }_{k=-\infty }^{+\infty }h(n)x(n-k)$

误差：$e(n)=d(n)-s^{\prime }(n)==d(n)-{\sum }_{k=-\infty }^{+\infty }h(n)x(n-k)$

均方误差：$J(h)=E[e^2(n)]$

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

$${\nabla }_{h}J=-2E[x(n-k)e(n)]$$

所以，
$$E[x(n-k)e(n)]=0(1)$$

$$E\left[x(n-k)\left(d(n)-\sum _{i=-\infty }^{\infty }h(i)x(n-i)\right)\right]=0$$

$$\sum _{i=-\infty }^{\infty }h(i)E\left[x(n-k)x(n-i)\right]=E\left[x(n-k)d(n)\right]$$

$$\sum _{i=-\infty }^{\infty }h(i)r_x(i-k)=r_{xd}(-k)$$

针对$M$阶FIR滤波器，维纳-霍夫方程（Wiener-Hopf）为${\sum }_{i=0}^{M-1}h(i)r_x(i-k)=r_{xd}(-k)$，那么，写成矩阵的形式就是
$$\begin{array}{l}Rh=r_{xd}\\ h=R^{-1}{r}_{xd}\end{array}$$
其中$R$是信号$x(n)$的自相关矩阵，$r_{xd}$是信号$x(n)和d(n)$的互相关矩阵。

仿真分析——Matlab代码

由于Wiener滤波器的参数求解过程中必须要用到参考信号，所以这里仿真采用的信号Signal为
$s=A\ast sin\left(2\pi f_{1}t\right)+B\ast sin\left(2\pi f_{2}t\right)$ ，
即为两个叠加的正弦信号。其中$A=10,B=15,f_1=1000,f_2=2000$ 。
根据采样定理，这里假定采样频率$f_s=100000$，采样间隔$T=1/f_s$，则$s_a(t)|{}_{t=nT}=s_a(nT)(0\le n\le 999)$。
对于不同的$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)$为原始信号$s(n)$
• 参考信号$d(n)$为加性高斯白噪声$v(n)$
• 参考信号$d(n)$为输入信号自身$x(n)$

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

一、参考信号$d(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: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(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)$为输入信号自身$x(n)$

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


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

总结

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