维纳滤波

维纳滤波 7.2 维 纳 滤 波 从连续的(或离散的)输入数据中滤除噪声和干扰以提取有用信息的过程称为滤波，而相应的装置称为滤波器。根据滤波器的输出是否为输入的线性函数，可将它分为线性滤波器和非线性滤波器两种。滤波器研究的一个基本课 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 就是:如何设计和制造最佳的或最优的滤波器。所谓最佳滤波器是指能够根据某一最佳准则进行滤波的滤波器。 20世纪40年代，维纳奠定了关于最佳滤波器研究的基础。即假定线性滤波器的输入为有用信号和噪声之和，两者均为广义平稳过程且知它们的二阶统计特性，维纳根据最小均方误差准则(滤波器的输出信号与需要信号之差的均方值最小)，求得了最佳线性滤波器的参数，这种滤波器被称为维纳滤波器。在维纳研究的基础上，人们还根据最大输出信噪比准则、统计检测准则以及其他最佳准则求得的最佳线性滤波器。实际上，在一定条件下，这些最佳滤波器与维纳滤波器是等价的。因而，讨论线性滤波器时，一般均以维纳滤波器作为参考。 维纳滤波理论用于解决最小均方误差下的线性滤波问题。设接收到(或观测到)的信号为随机信号 (7-1) 其中s(t)是未知的实随机信号，n(t)是噪声。要设计的线性滤波器，其冲击响应为h(t, τ)，输入为x(t)，输出为，即 (7-2) 令为估计误差。冲击响应h(t, τ)按最小均方误差准则确定，即h(t, τ)必须满足使 (7-3) 达到最小。根据最小均方误差估计的正交条件，有以下关系成立 (7-4) 令 (7-5) (7-6) 则有 (7-7) 上述方程通常称为非平稳随机过程条件下的维纳-霍甫(Wiener-Kolmogorov)积分方程。特别当x(t)，s(t)均为广义(或宽)平稳随机信号，而滤波器是线性时不变系统的情况下，x(t)与s(t)必为联合平稳，式(7-7)可写为 (7-8) 令，，则有 (7-9) 此处，“*”号表示卷积，对上式两边取Fourier变换，可得 (7-10) (7-11) 对于因果线性系统，有 (7-12) 采用完全相同的分析方法，推得因果平稳维纳-霍甫积分方程如下 (7-13) (7-14) 其中，表示的零、极点位于， 表示的零、极点位于。表示位于的零、极点。 MATLAB图像处理工具箱提供了wiener2函数进行自适应滤出图像噪声，它根据图像的局 部方差来调节滤波器的输出，往往较线性滤波效果好，可以更好地保存图像的边缘和高频细节信息。 Wiener2函数采用的算法是首先估计像素的局部均值和方差: (7-15) Ø (7-16) 其中，Ω是图像中每个像素的M×N的邻域。然后，对每个像素利用wiener滤波器估计其灰度值: (7-17) 2式中，v是图像中噪声的方差。 Wiener2的语法格式为: J=wiener2(I,[m,n]) J=wiener2(I,[m,n],noise) [J,noise]=wiener2(I,[m,n]) 其中，J=wiener2(I,[m,n])返回有噪声图像I经过wierner(维纳)滤波后的图像，[m,n]指定滤波器窗口大小为，默认值为。J=wiener2(I,[m,n],noise)指定噪声的功率，[J,noise]=wiener2(I,[m,n])在图像滤波的同时，返回噪声功率的估计值noise。 【例7-1】对加入高斯噪声的图像saturan.png作维纳滤波。 例程7-1 噪声图像维纳滤波 % e.g.7-1.m for example7-1; %test the function of weina filter. RGB = imread('saturn.png'); I = rgb2gray(RGB); J = imnoise(I,'gaussian',0,0.005); figure, imshow(J); K = wiener2(J,[5 5]); %指定滤波器窗口大小 figure, imshow(K); 图7-1 噪声图像 图7-2 维纳滤波复原图像 实现维纳滤波的要求是:?输入过程是广义平稳的;?输入过程的统计特性是已知的。根据其他最佳准则的滤波器亦有同样要求。然而，由于输入过程取决于外界的信号、干扰环境，这种环境的统计特性常常是未知的、变化的，因而难以满足上述两个要求。这就促使人们研究自适应滤波器。 卡尔曼滤波器 由于维纳滤波器难以用于实时处理，满足不了20世纪50年代航天航空技术发展的要求，于是人们开始探索新的理论和技术途径。20世纪60年代新出现的卡尔曼滤波理论，用信号与噪声的状态空间模型取代了相关函数，用时域的微分方程来表示滤波问题，得到了递推估计算法，适用于计算机实时处理，它突破了维纳滤波的平稳过程的限制，也没有无限时间的要求，这一对维纳滤波理论的重大突破很快地被用于空间技术、自动控制和信号处理等领域。卡尔曼滤波由滤波方程和预测方程两部分组成。 7.5.1 滤波基本方程 设信号状态方程和量测方程分别为 (7-133) (7-134) 式中，为信号的状态向量，为量测向量，和分别为状态噪声和量测噪声，且为互不相关的高斯白噪声向量序列，其协方差分别为和 ;，和分别为状态转移矩阵、输入矩阵和观测矩阵。 卡尔曼滤波基本方程为 (7-135) (7-136) (7-137) (7-138) (7-139) 其中，残差(新息)定义为 (7-140) 协方差矩阵为 (7-141) 7.5.2 一步预测基本方程 卡尔曼一步预测基本方程为 (7-142) (7-143) (7-144) 式中，为一步预测增益阵。 【例7-16】利用卡尔曼滤波器跟踪机动目标。 例程7-14 基于卡尔曼滤波器的机动目标跟踪 % Make a point move in the 2D plane % State = (x y xdot ydot). We only observe (x y). % X(t+1) = Φ(t) X(t) + noise(Q) % Y(t) = H X(t) + noise(R) ss = 4; % state size os = 2; % observation size F = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1]; H = [1 0 0 0; 0 1 0 0]; Q = 0.1*eye(ss); R = 1*eye(os); initx = [10 10 1 0]'; %target initial parameters initV = 10*eye(ss); seed = 9; rand('state', seed); randn('state', seed); T = 15; [x,y] = sample_lds(F, H, Q, R, initx, T); %generate target data %kalman filter [xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y, F, H, Q, R, initx, init V); % one step predict [xsmooth, Vsmooth] = kalman_smoother(y, F, H, Q, R, initx, initV); %calculate the error between the filtered data and the real data dfilt = x([1 2],:) - xfilt([1 2],:); mse_filt = sqrt(sum(sum(dfilt.^2))); dsmooth = x([1 2],:) - xsmooth([1 2],:); mse_smooth = sqrt(sum(sum(dsmooth.^2))) figure(1) clf %subplot(2,1,1) hold on plot(x(1,:), x(2,:), 'ks-'); plot(y(1,:), y(2,:), 'g*'); plot(xfilt(1,:), xfilt(2,:), 'rx:'); for t=1:T, plotgauss2d(xfilt(1:2,t), Vfilt(1:2, 1:2, t)); end hold off legend('true', 'observed', 'filtered', 3) xlabel('x') ylabel('y') % 3x3 inches set(gcf,'units','inches'); set(gcf,'PaperPosition',[0 0 3 3]) %print(gcf,'-depsc','/home/eecs/murphyk/public_html/Bayes/Figures/aima_filter ed.eps'); %print(gcf,'-djpeg','-r100', '/home/eecs/murphyk/public_html/Bayes/Figures/ai ma_filtered.jpg'); figure(2) %subplot(2,1,2) hold on plot(x(1,:), x(2,:), 'ks-'); plot(y(1,:), y(2,:), 'g*'); plot(xsmooth(1,:), xsmooth(2,:), 'rx:'); for t=1:T, plotgauss2d(xsmooth(1:2,t), Vsmooth(1:2, 1:2, t)); end hold off legend('true', 'observed', 'smoothed', 3) xlabel('x') ylabel('y') % 3x3 inches set(gcf,'units','inches'); set(gcf,'PaperPosition',[0 0 3 3]) %print(gcf,'-djpeg','-r100', '/home/eecs/murphyk/public_html/Bayes/Figures/ai ma_smoothed.jpg'); %print(gcf,'-depsc','/home/eecs/murphyk/public_html/Bayes/Figures/aima_smoo thed.eps'); 图7-27 卡尔曼滤波器跟踪机动目标仿真 图7-27中“?”表示机动目标的真实轨迹，“,”表示目标观测轨迹，“…”表示经卡尔曼滤波后的 目标轨迹，圆表示跟踪波门。可以看出，卡尔曼滤波器能较为准确地跟踪机动目标。函数k alman_filter和 kalman_update分别是实现卡尔曼滤波方程滤波和状态更新的代码。 function [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, varar gin) % Kalman filter. % [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, ...) % % INPUTS: % y(:,t) - the observation at time t % A - the system matrix % C - the observation matrix % Q - the system covariance % R - the observation covariance % init_x - the initial state (column) vector % init_V - the initial state covariance % % OPTIONAL INPUTS (string/value pairs [default in brackets]) % 'model' - model(t)=m means use params from model m at time t [ones (1,T) ] % In this case, all the above matrices take an additional final dimension, % i.e., A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m). % However, init_x and init_V are independent of model(1). % 'u'- u(:,t) the control signal at time t [ [] ] % 'B'- B(:,:,m) the input regression matrix for model m % % OUTPUTS (where X is the hidden state being estimated) % x(:,t) = E[X(:,t) | y(:,1:t)] % V(:,:,t) = Cov[X(:,t) | y(:,1:t)] % VV(:,:,t) = Cov[X(:,t), X(:,t-1) | y(:,1:t)] t >= 2 % loglik = sum{t=1}^T log P(y(:,t)) % % If an input signal is specified, we also condition on it: % e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t)] % If a model sequence is specified, we also condition on it: % e.g., x(:,t) = E[X(:,t) | y(:,1:t), u(:, 1:t), m(1:t)] [os T] = size(y); ss = size(A,1); % size of state space % set default params model = ones(1,T); u = []; B = []; ndx = []; args = varargin; nargs = length(args); for i=1:2:nargs switch args{i} case 'model', model = args{i+1}; case 'u', u = args{i+1}; case 'B', B = args{i+1}; case 'ndx', ndx = args{i+1}; otherwise, error(['unrecognized argument ' args{i}]) end end x = zeros(ss, T); V = zeros(ss, ss, T); VV = zeros(ss, ss, T); loglik = 0; for t=1:T m = model(t); if t==1 %prevx = init_x(:,m); %prevV = init_V(:,:,m); prevx = init_x; prevV = init_V; initial = 1; else prevx = x(:,t-1); prevV = V(:,:,t-1); initial = 0; end if isempty(u) [x(:,t), V(:,:,t), LL, VV(:,:,t)] = ... kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, 'initial', initial); else if isempty(ndx) [x(:,t), V(:,:,t), LL, VV(:,:,t)] = ... kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, ... 'initial', initial, 'u', u(:,t), 'B', B(:,:,m)); else i = ndx{t}; % copy over all elements; only some will get updated x(:,t) = prevx; prevP = inv(prevV); prevPsmall = prevP(i,i); prevVsmall = inv(prevPsmall); [x(i,t), smallV, LL, VV(i,i,t)] = ... kalman_update(A(i,i,m), C(:,i,m), Q(i,i,m), R(:,:,m), y(:,t), prevx(i), prev Vsmall, ... 'initial', initial, 'u', u(:,t), 'B', B(i,:,m)); smallP = inv(smallV); prevP(i,i) = smallP; V(:,:,t) = inv(prevP); end end loglik = loglik + LL; end function [xnew, Vnew, loglik, VVnew] = kalman_update(A, C, Q, R, y, x, V, varargin) % KALMAN_UPDATE Do a one step update of the Kalman filter % [xnew, Vnew, loglik] = kalman_update(A, C, Q, R, y, x, V, ...) % % INPUTS: % A - the system matrix % C - the observation matrix % Q - the system covariance % R - the observation covariance % y(:)- the observation at time t % x(:) - E[X | y(:, 1:t-1)] prior mean % V(:,:) - Cov[X | y(:, 1:t-1)] prior covariance % % OPTIONAL INPUTS (string/value pairs [default in brackets]) % 'initial' - 1 means x and V are taken as initial conditions (so A and Q ar e ignored) [0] % 'u'- u(:) the control signal at time t [ [] ] % 'B'- the input regression matrix % % OUTPUTS (where X is the hidden state being estimated) % xnew(:) =E[ X | y(:, 1:t) ] % Vnew(:,:) = Var[ X(t) | y(:, 1:t) ] % VVnew(:,:) = Cov[ X(t), X(t-1) | y(:, 1:t) ] % loglik = log P(y(:,t) | y(:,1:t-1)) log-likelihood of innovatio % set default params u = []; B = []; initial = 0; args = varargin; for i=1:2:length(args) switch args{i} case 'u', u = args{i+1}; case 'B', B = args{i+1}; case 'initial', initial = args{i+1}; otherwise, error(['unrecognized argument ' args{i}]) end end % xpred(:) = E[X_t+1 | y(:, 1:t)] % Vpred(:,:) = Cov[X_t+1 | y(:, 1:t)] if initial if isempty(u) xpred = x; else xpred = x + B*u; end Vpred = V; else if isempty(u) xpred = A*x; else xpred = A*x + B*u; end Vpred = A*V*A' + Q; end e = y - C*xpred; % error (innovation) n = length(e); ss = length(A); S = C*Vpred*C' + R; Sinv = inv(S); ss = length(V); loglik = gaussian_prob(e, zeros(1,length(e)), S, 1); K = Vpred*C'*Sinv; % Kalman gain matrix % If there is no observation vector, set K = zeros(ss). xnew = xpred + K*e; Vnew = (eye(ss) - K*C)*Vpred; VVnew = (eye(ss) - K*C)*A*V;