Wiener維納濾波基本原理及其演算法實現
文章轉載自：http://blog.sina.com.cn/s/blog_bb81c2230102xdbl.html 如果有侵權，請聯絡博主刪除
To learn, to share, to debate, then comes progress.
1.演算法背景：
訊號濾波的實質為從觀測訊號中提取有效訊號，隨著數學理論的發展與實際應用的需求，基於不同原理的濾波方法被不斷地提出來，雖然依據的準則，推導的過程各有差異，但最終的目的均是減小訊號估計的誤差，使濾波系統的輸出訊號儘可能地接近實際訊號。
Wiener濾波是第二次世界大戰中，為了解決火力控制系統精確跟蹤問題，Wiener相繼提出了平穩隨機過程的最優線性濾波理論，首次將數理統計知識和線性系統理論聯絡起來，形成了對隨機訊號作平滑，濾波和預測的最新估計理論。在此後的發展中，Wiener濾波被應用於更多的領域，並沿用至今。
2.演算法原理：
（1）有限長濾波器
對於一列輸入訊號x，一般的無限長線性濾波器輸出為：
y(n)= Σh(m)x(n-m) m=0…∞
實際中，濾波器的長度，即階數是有限長的，設為M，則有：
y(n)= Σh(m)x(n-m) m=0…M
即濾波器的當前時刻輸出為前M個時刻的值經過加權之後得到的。
為便於書寫與理解，上式可以寫為矩陣形式：
y(n)=H(m)*X(n)
如果期望訊號d已知，則可以計算輸出與期望訊號之間的誤差：
e(n)=d(n)-y(n)= d(n)- H(m)*X(n) m=0…M
Wiener濾波的目標就是，如何確定一個長為M的係數序列H，使得上述誤差值最小。
（2）最小均方誤差濾波
根據目標函式的不同，又可以將濾波演算法細分為不同的類別，一般來說有最小均方誤差，最小二乘誤差等等，這裡只討論最小均方誤差。
令目標函式為：
Min E[e(n)^2]= E[(d(n)- H(m)*X(n))^2]
當濾波器的係數最優時，目標函式對係數的倒數應該為0，即：
dE[e(n)^2]/dH=0
2 E[ (d(n)- H(m)*X(n))]* X(n)=0
E[(d(n) X(n))- H(m)E[X(n)X(n)]=0
根據隨機過程的知識，上式可以表達為：
Rxd-H*Rxx=0
其中Rxd與Rxx分別為輸入訊號與期望訊號的相關矩陣與輸入訊號的自相關矩陣。
從而有：
H=Rxx-1*Rxd
至此，便得到了Wiener濾波的基本原理與公式推導。
3.演算法應用與實現
理解了演算法的原理之後，下邊舉一個小的例子來考察如何應用Wienar濾波處理實際問題。
問題背景：一個點目標在x，y平面上繞單位圓做圓周運動，由於外界干擾，其運動軌跡發生了偏移。其中，x方向的干擾為均值為0，方差為0.05的高斯噪聲；y方向干擾為均值為0，方差為0.06的高斯噪聲。
問題分析與思路：
將物體的運動軌跡分解為X方向和Y方向，並假設兩個方向上運動相互獨立。分別將運動軌跡離散為一系列點，作為濾波器的輸入，分別在兩個方向上進行濾波，最終再合成運動軌跡。
程式設計思路：
生成期望訊號-新增噪聲-計算相關矩陣-求解最佳濾波器係數-濾波運算-輸出訊號-合成軌跡
4.結果與分析
5.原始碼
%***********************************************
%該程式使用Wiener濾波方法對圓周運動軌跡進行控制
%訊號模型：d=s+no 觀測訊號=期望訊號+噪聲訊號
%進行一次Wiener濾波，得到最佳濾波器係數
17.4 by Howie
clear
close all
N=500;
theta=linspace(0,2*pi,N); %極座標引數
s_x=cos(theta); %x，y方向上的期望訊號
s_y=sin(theta);
no_x=normrnd(0,sqrt(0.05),1,N); %高斯白噪聲
no_y=normrnd(0,sqrt(0.06),1,N);
d_x=s_x+no_x; %觀測訊號
d_y=s_y+no_y;
M=500;%M為濾波器的階數
%% 對x方向上資料進行濾波
rxx=xcorr(d_x);
Rxx=zeros(N);
% temp=toeplitz(rxx);
for i=1:N %觀測訊號的相關矩陣
for j=1:N
Rxx(i,j)=rxx(N+i-j);
end
end
rxd=xcorr(s_x,d_x); %觀測訊號與期望訊號的相關矩陣
Rxd=rxd(N:N+M-1); %向量而非矩陣
hopt_x=Rxx\Rxd';
% de_x=conv(hopt_x,d_x);
de_x=zeros(1,N);
for n=1:N
for i=1:n-1
de_x(n)=de_x(n)+hopt_x(i)*d_x(n-i);
end
end
de_x(1:2)=d_x(1:2);
ems_x=sum(d_x.^2)-Rxd*hopt_x;
e_x=de_x-s_x;
% de_x(N-1:N)=d_x(N-1:N);
%% 對y方向上資料進行濾波 處理思路同x方向
ryy=xcorr(d_y);
Ryy=zeros(N);
for i=1:N
for j=1:N
Ryy(i,j)=ryy(N+i-j);
end
end
% temp=toeplitz(ryy);
% Ryy=temp(1:M,N:N+M-1);
ryd=xcorr(s_y,d_y);
% temp=toeplitz(ryd);
% Ryd=temp(1:N,N:length(temp));
Ryd=ryd(N:N+M-1);
hopt_y=Ryy\Ryd';
% de_y=conv(hopt_y,d_y);
de_y=zeros(1,N);
for n=1:N
for i=1:n-1
de_y(n)=de_y(n)+hopt_y(i)*d_y(n-i);
end
end
de_y(1:2)=d_y(1:2);
ems_y=sum(d_y.^2)-Ryd*hopt_y;
e_y=de_y-s_y;
% de_y(N-1:N)=d_y(N-1:N);
%% plot
figure
plot(s_x,s_y,'r','linewidth',2)
hold on
plot(d_x,d_y,'b')
hold on
plot(de_x,de_y,'k-')
title('維納濾波預測軌跡')
legend('期望軌跡','觀測軌跡','濾波軌跡')
%% %% x方向上繪圖
figure
suptitle('X方向上維納濾波效果')
subplot(321)
plot(s_x)
title('期望訊號')
subplot(322)
plot(no_x)
title('噪聲訊號')
subplot(323)
plot(d_x)
title('觀測訊號')
subplot(324)
plot(de_x)
title('濾波後訊號')
subplot(325)
plot(ems_x,'o')
title('最小均方誤差')
subplot(326)
plot(e_x)
title('絕對誤差')
%% y方向上繪圖
figure
suptitle('Y方向上維納濾波效果')
subplot(321)
plot(s_y)
title('期望訊號')
subplot(322)
plot(no_y)
title('噪聲訊號')
subplot(323)
plot(d_y)
title('觀測訊號')
subplot(324)
plot(de_y)
title('濾波後訊號')
subplot(325)
plot(ems_y,'o')
title('最小均方誤差')
subplot(326)
plot(e_y)
title('絕對誤差')