维纳滤波及其简单实现

文章目录

介绍基本概念简单实现过程matlab实现结果小结

介绍

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

基本概念

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

自相关函数 为了描述随机变量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)$互相关函数 描述不同随机过程引入的随机序列的相互关系 $$R_{XY}(t)=E[X(n)\ast Y(n+t)]$$ 并且满足：$R_{XY}(t)=R_{YX}(-t)$维纳滤波原理 本质上，维纳滤波是设计一个冲击响应函数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矩阵的时候利用了这一点。

简单实现过程

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

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函数计算得到，这里为了完整地了解整个过程计算过程，使用自己编写的代码。