维纳滤波及其简单实现

介绍

随机信号包括了确定信号和随机噪声两部分。维纳滤波的本质是设计一组冲击响应的函数，抑制信号中的随机噪声部分，或者说非预期信号部分，使得信号与预期值的均方误差达到最小。

基本概念

在开始维纳滤波的介绍前，先描述一下几个基本的概念
以下只给出离散过程的公式

1. 自相关函数
   为了描述随机变量X(n),X(n+t)，在不同时刻下的相互联系，引入了自相关函数
   t为间隔,t>0。
   $$R_{xx}(t)=E[X(n)\ast X(n+t)]$$
   作为一种预期，这里自相关函数只与间隔有关，与起始时间无关。同时，要注意到，自相关函数是关于时间间隔t的偶函数，既有$R_{xx}(t)=R_{xx}(-t)$
2. 互相关函数
   描述不同随机过程引入的随机序列的相互关系
   $$R_{XY}(t)=E[X(n)\ast Y(n+t)]$$
   并且满足：$R_{XY}(t)=R_{YX}(-t)$
3. 维纳滤波原理
   本质上，维纳滤波是设计一个冲击响应函数h，对观察到的信号序列X进行滤波，使得其与期望信号S的最小平方和达到最小值。
   假设信号S的长度为N，下标从0开始计算,$S=[s(0),\cdots ,s(N-1)],s(i)\in \mathbb{R}$,
   对于滤波器h而言，其滤波过程如下，
   对n时刻的信号估计值
   $$\hat{s}(n)=\sum _{m=0}^{n}h(m)x(n-m)$$
   维纳滤波器的目标如下：
   $$minE(e_n^2)=E([s(n)-\hat{s}(n){]}^2)$$
   按$h(i)$对上式进行求导，由最小值条件得到
   $$\begin{gathered} E([s(n)-\hat{s}(n)]\ast x(n-i))=0\Rightarrow \\ E[s(n)x(n-i)]-E[\sum _{m=0}^{n}h(m)x(n-m)\ast x(n-i)]=0 \end{gathered}$$
   整理上式，结合相关函数的性质可以得到
   $$\begin{gathered} E[s(n)x(n-i)]=R_{sx}(-i)=R_{xs}(i) \\ E[h(m)x(n-m)\ast x(n-i)]=h(m)E(x(n-m)\ast x(n-i))=h(m)R_{xx}(m-i) \\ {R}_{xs}(i)=\sum _{m=0}^{n}h(m)R_{xx}(m-i) \end{gathered}$$
   对于长度为N的滤波器h，则有
   $$\left[\begin{array}{c}{R}_{xs}(0)\\ ⋮\\ R_{xs}(i)\\ ⋮\\ R_{xs}(N-1)\end{array}\right]=\left[\begin{array}{ccccc}{R}_{xx}(0)& \cdots & R_{xx}(i)& \cdots & R_{xx}(N-1)\\ ⋮& \ddots & ⋮& & ⋮\\ R_{xx}(-i)& \cdots & R_{xx}(0)& \cdots & R_{xx}(N-1-i)\\ ⋮& & ⋮& \ddots & ⋮\\ R_{xx}(1-N)& \cdots & R_{xx}(i-N+1)& \cdots & R_{xx}(0)\end{array}\right]\left[\begin{array}{c}h(0)\\ ⋮\\ h(i)\\ ⋮\\ h(N-1)\end{array}\right]$$
   将上式进行简化，得到 $RXS=RXX\ast h\Rightarrow h=RX{X}^{-1}RXS$
   $RXX$是一个对称矩阵，且对角线元素都是$R_{xx}(0)$,在下面构造RXX矩阵的时候利用了这一点。

简单实现过程

1. 构造正弦信号s
2. 加入白噪声noise
3. 合成最终信号 x = s+noise
4. 生成RXX矩阵和RXS向量
5. 计算h滤波器
6. 信号还原并显示计算结果

matlab实现

clear;
N=600;               %data size and filter size
theta=linspace(0,2*pi,N);
s=sin(theta);
noise=normrnd(0,sqrt(0.05),1,N);    %noise
x = s+noise;
RXX = ConstructRxxMatrix(x);
RXS = ConstructRxsVector(x,s);
h = RXX\RXS;
%recovery s1 from x using h filter
s1 = zeros(N,1);
for i = 1:1:N
    frag1 = h(1:i);frag2 = x(i:-1:1);
    s1(i) = dot(frag1,frag2);
end
subplot(2,2,1);plot(s);title('expected signal');
subplot(2,2,2);plot(x);title('real signal');
subplot(2,2,3);plot(noise);title('white noise');
subplot(2,2,4);plot(s1);title('denoise signal');
%Calc Rxs(t)  && only handle t>0 
function r = RelateValue(x,s,t)
n  = length(x);
frag1 = s(1+t:n);frag2 = x(1:(n-t));
r = dot(frag1,frag2);
end
%Construct Rxx Matrix%
function Rxx = ConstructRxxMatrix(x)
n = length(x);
Rxx = zeros(n);
RV = zeros(n,1);
for i = 1:1:n
    RV(i)  = RelateValue(x,x,i-1);
end
for i = 1:1:n
    for j = i+1:1:n
        Rxx(i,j) = RV(j-i+1);
    end
end
Rxx = Rxx+Rxx'+diag(ones(n,1))*RV(1);
end
function Rxs = ConstructRxsVector(x,s)
n = length(x);
Rxs = zeros(n,1);
for i = 1:1:n
    Rxs(i) = RelateValue(x,s,i-1);
end
end

结果

小结

和卡尔曼滤波一样，维纳滤波是信号处理中一种经典的滤波算法。上述的互相关函数可以借用matlab的xcorr函数计算得到，这里为了完整地了解整个过程计算过程，使用自己编写的代码。